Abstract:  Icejams  cause  flooding  in  northern  temperate- 
climate  areas,  usually  forming  rapidly,  often  with  little 
warning,  constricting  water  flow  and  elevating  water 
levels.  Consequently,  jam  formation  comprises  highly 
unsteady  processes:  drifting  ice  pieces  are  brought  to 
rest,  accumulated  ice  shoves  and  thickens,  and  initial 
water  depths  and  velocities  change.  Those  processes 
are  even  more  unsteady  when  a  jam  collapses.  Prior 
simulations  of  icejams,  however,  treat  them  as  simply 
stationary,  uniformly  thick  accumulations  of  ice  pieces. 
No  account  is  taken  of  the  impact  forces  exerted  by 
moving  ice,  an  estimation  that  is  further  complicated  by 
the  need  to  couple  equations  describing  water  flow  and 
ice  movement.  Under  the  dynamic  conditions  attendant 


tojam  formation,  waterflowand  ice  movement  interactively 
influence  each  other.  This  report  evaluates  the  importance 
of  ice  momentum  on  ice  jam  thickness  and  thickness 
distribution  using  experiments  conducted  with  laboratory 
flumes  and  a  numerical  model  in  which  the  equations  of 
motion  for  one-dimensional  flow  of  water  and  ice  are 
solved  as  fully  coupled.  In  this  regard,  the  model  is 
unique,  enabling  simulation  of  the  important  unsteady 
interactions  of  water  and  ice,  and  determination  of  their 
effects  on  jam  thickness.  Ice  momentum  should  be  taken 
into  accountfor  mostjams  because  it  leads  to  significantly 
thicker  jams  and  affects  the  thickness  profile.  A  useful 
dimensionless  parameter  is  identified  for  generalizing  this 
finding. 
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Unsteady  Ice  Jam  Processes 

JON  E.  ZUFELT  AND  ROBERT  ETTEMA 


INTRODUCTION 

Background 

Ice  jams  cause  massive  damage  annually  throughout  the  world's  northern  tem¬ 
perate  regions.  In  the  U.S.  alone,  the  annual  damages  from  ice  jams  average  $125 
million  (USACE  1994).  These  include  property  losses,  emergency  assistance,  flood 
insurance,  and  increased  operation  and  maintenance  costs,  replacement  of  infra¬ 
structure,  and  loss  of  hydropower  revenues.  Most  damages  are  caused  by  high 
water  levels  associated  with  ice  jams,  though  some  are  from  the  direct  impact  of 
moving  ice  during  ice  runs. 

The  formation  and  evolution  of  ice  jams  comprise  a  series  of  inherently  unsteady 
processes,  in  which  moving  ice  is  brought  to  rest  in  accumulations  that  shove  and 
thicken  in  accordance  with  changing  forces  exerted  by  water  flow,  accumulation 
weight,  and  bank  roughness.  These  processes  are  even  more  unsteady  when  a  jam 
collapses,  plows  downstream,  and  possibly  reforms.  Prior  formulations  for  ice  jams 
treat  them  as  stationary,  steady-state  ice  accumulations  that  are  subject  to  invariant 
flow  conditions.  This  study,  however,  presents  the  first  formulation  for  and  exami¬ 
nation  of  the  fully  coupled  dynamic  nature  of  the  unsteady  processes  associated 
with  jam  formation. 

Need  for  research 

In  efforts  to  protect  life  and  property  from  the  damages  attendant  to  ice  jams 
and  related  flooding,  prior  models  were  developed  to  predict  water-level  changes 
caused  by  ice  jams.  Those  models  treat  the  evolution  of  ice  jam  thickness  (shoving 
and  thickening)  as  quasi-steady,  with  jam  thickness  spontaneously  adjusting  to  a 
new  equilibrium  value  in  concert  with  water  flow  changes.  Steady-flow  models, 
such  as  HEC-2  modified  with  ice  cover  option,  simply  provide  the  steady  water 
levels  that  would  exist  with  a  uniformly  thick  jam  already  in  place.  The  long-stand¬ 
ing  assumption  used  is  of  an  ice  jam  of  equilibrium  thickness,  floating  in  static 
force  equilibrium  just  on  the  verge  of  stability  or  failure.  Other  models  simulate  the 
unsteadiness  of  the  water  flow  using  the  conservation  of  mass  and  momentum 
equations  for  the  water,  but  solve  in  an  uncoupled  manner  for  thickness  between 
time  steps,  again  by  the  static  force  balance. 

There  currently  are  no  formulations  that  describe  the  coupled  interaction  of  the 
water  and  ice  movement  and  their  effects  on  flow  depth  and  ice  thickness.  Also,  no 
information  exists  on  how  jams  evolve,  fail,  and  thicken.  In  that  regard,  the  follow¬ 
ing  important  groups  of  questions  need  to  be  addressed: 

•  How  do  ice  jams  evolve?  Present  formulations  allow  for  instantaneous  changes 
in  the  jam  thickness  attributable  to  changes  in  the  forces  acting  on  the  jam.  No 
account  is  currently  made  for  the  impact  forces  generated  by  moving  ice.  Once 
a  jam  fails,  ice  is  mobilized  and  travels  downstream,  often  at  high  speed 
(Henderson  and  Gerard  1981).  Do  jams  move  and  then  thicken  upon  failure, 
thicken  as  they  fail,  or  thicken  and  result  in  a  progressive  downstream-mov¬ 
ing  failure? 


•  To  what  extent  does  ice  momentum  affect  jam  thickness?  If  jams  move  upon 
failure,  the  force  levels  acting  on  the  jam  inevitably  change.  What  are  these 
changes  and  how  do  they  affect  the  forces  acting  on  the  jam? 

•  What  are  the  effects  of  the  interaction  of  the  water  and  ice  motion?  The  water 
shear  stress  on  the  jam  underside  is  one  of  the  principal  forces  on  the  jam. 
When  a  jam  fails  and  moves,  however,  this  force  is  reduced,  which  interac¬ 
tively  reduces  the  resistance  to  water  flow.  As  water  and  ice  motion  are  inter¬ 
related,  what  are  the  consequences,  for  jam  thickness  prediction,  of  uncou¬ 
pling  their  influences  as  is  currently  done  in  existing  formulations? 

Field  observations 

The  inherent  unsteadiness  of  jam  formation  and  failure  is  obvious  from  field 
observations.  As  an  example,  observations  of  a  freezeup  jam  on  the  Salmon  River 
near  Salmon,  Idaho,  made  from  a  small  bridge  about  1  km  downstream  of  the  lead¬ 
ing  edge  (head)  of  the  jam,  are  presented  here.  The  jam  was  approximately  25  km 
long,  and,  owing  to  mild  weather  (-3°C),  the  leading  edge  had  been  stationary 
over  the  night.  Water  levels  were  slightly  more  than  bankfull  at  the  bridge  and 
thick  frazil  accumulations  filled  the  channel. 

The  water  level  began  rising,  first  noticeably  at  the  treeline  along  the  bank,  then 
in  mid-channel,  as  the  water  began  to  flood  the  surface  of  the  jam  (Fig.  1).  The  ice  in 
the  channel  appeared  to  rise  slightly,  and  shear  cracks  could  be  seen  forming  about 
10  m  out  from  the  banks.  The  water  levels  continued  to  rise  and  the  ice,  groaning, 
began  to  slowly  move  downstream  en  masse  (Fig.  2).  In  a  matter  of  seconds,  the 
entire  channel  section  of  the  jam  for  about  1  km  on  either  side  of  the  bridge 
(between  the  shear  cracks)  was  moving  downstream  (Fig.  3).  As  the  ice  moved,  the 
water  level  fell,  until  the  center  portion  of  the  channel  was  clear  of  ice  upstream  to 
where  the  leading  edge  had  previously  been  located.  Once  the  ice  had  passed, 
water  levels  dropped  by  approximately  0.5  m,  exposing  shear  walls  of  ice  along  the 


Figure  2.  Ice  jam  failure  at  beginning  of  ice  motion;  view  is  looking  upstream. 


banks  (Fig.  4).  The  shear  walls  were  grounded  on  the  bed.  The  unsteadiness 
observed  is  significant  because  the  time  lapsed  since  the  water  levels  initially  rose 
to  the  final  passage  of  the  ice  from  upstream  was  only  about  15  minutes.  Questions 
remain  as  to  what  happened  to  the  ice  as  it  traveled  downstream,  what  caused  the 
initial  water  level  rise,  and  what  combination  of  water  and  ice  flow  resulted  in  the 
initial  accumulation. 


Figure  3.  Ice  jam  failure  with  ice  fully  mobilized;  view  is  looking  downstream. 
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Figure  4.  River  channel  following  ice  jam  failure;  view  is  looking  upstream. 


Objective  and  approach 

The  ultimate  objective  of  this  work  is  to  address  the  three  sets  of  questions  pre¬ 
sented  in  the  Needs  for  Research  section.  To  do  so  required  carrying  out  the  follow¬ 
ing  tasks: 

•  Determine  the  temporal  sequence  of  events  associated  with  ice  jam  formation, 
failure,  and  reformation. 

•  Identify  the  important  processes  and  parameters  involved  in  shoving  and  thick¬ 
ening,  and  properly  formulate  the  equations  describing  them. 

•  Develop  a  numerical  model,  correctly  representing  the  shoving  and  thicken¬ 
ing  of  ice,  for  use  in  examining  the  formation  and  evolution  of  jams,  including 
freezeup  and  breakup  jams.  The  model  would  be  used  further  to  investigate 
1)  the  progression  of  freezeup  jams  formed  at  an  ice  boom;  2)  the  effects  of  a 
jam  on  hydraulics  as  the  length  of  the  jam  increases  and  evolves  by  shoving 
and  thickening;  3)  the  failure  of  an  ice  jam  ascribable  to  increases  in  water 
discharge  simulating  the  effects  of  hydropower  releases  or  surges  from  the 
failure  of  upstream  ice  jams;  and  4)  the  effects  of  ice  momentum  on  the  pre¬ 
dicted  thickness  of  jams. 

•  Identify  a  parameter  that  delineates  when  the  effects  of  ice  momentum  be¬ 
come  important  for  determining  a  jam  thickness  profile,  and  when  a  fully 
coupled,  moving  ice  model  should  be  used  instead  of  a  steady-state,  station¬ 
ary  ice  model. 

Discussion 

Although  the  published  literature  on  jams  contains  descriptions  of  the  general 
processes  leading  to  shoving  and  thickening,  there  is  no  description  of  how  jams  in 
fact  move  and,  in  so  doing,  modify  water  flow.  In  very  general  terms,  an  ice  jam 
forms  when  the  downstream  movement  of  ice  is  stopped.  If  the  forces  exerted  on 
the  jam  continue  to  increase  (owing  to  increased  water  discharge,  increased  cover 


length,  or  reduced  jam  strength),  the  jam  eventually  fails  and  moves  downstream. 
If  an  area  of  channel  downstream  is  encountered  where  the  resisting  forces  on  the 
moving  ice  are  again  great  enough,  the  jam  will  reform.  Existing  models  predict 
the  equilibrium  jam  thickness,  which  is  the  constant  thickness  that  would  be 
expected  in  a  uniform  channel  under  conditions  of  steady  flow  when  the  resisting 
and  downstream-acting  forces  are  perfectly  in  balance.  Those  models  assume  that, 
when  the  net  downstream-acting  force  reaches  the  level  of  passive  pressure  failure, 
the  jam  must  thicken  to  withstand  the  forces.  The  unsteadiness  of  both  the  ice  and 
water  movement  during  a  shoving  and  thickening  event  make  the  concept  of  equi¬ 
librium  thickness  questionable. 

In  contrast  with  jam  formation,  juxtaposition  (or  surface  assembly)  of  ice  floes  is 
primarily  a  single-layer  process  that  can  be  adequately  described  from  hydraulics 
considerations,  and  from  the  size,  shape,  and  distribution  of  ice  pieces.  Much  work 
has  addressed  the  problem  of  block  undertuming  at  the  upstream  edge  of  obstruc¬ 
tions,  and  some  work  has  addressed  what  happens  to  the  blocks  following 
undertuming.  Such  cover-formation  processes  are  easy  to  visualize  in  the  labora¬ 
tory  and  the  river,  as  they  generally  occur  at  the  water  surface  and  entail  the 
motion  of  single-layer  ice  floes  coming  into  contact  with  a  stationary  obstruction. 
The  juxtaposition  of  ice  floes  or  the  motion  of  an  ice  block  at  the  upstream  edge  of 
an  obstruction  can  be  considered  to  be  a  fairly  steady  process  because  the  effects  of 
ice-piece  movement  on  the  hydraulics  are  minimal. 

Shoving  and  thickening,  however,  are  much  more  common  in  nature  during  the 
development  of  freezeup  jams  (made  up  of  frazil  slush  or  pans  and  small  ice  pieces), 
as  well  as  during  the  formation  and  evolution  of  breakup  jams.  The  manner  whereby 
jams  form  and  evolve  is  also  important  in  determining  how  they  fail.  In  compari¬ 
son  to  the  fairly  steady  water  and  ice  motion  during  juxtaposition  or  undertuming, 
the  ice  and  water  interaction  during  jam  failure  and  thickening  results  in  highly 
unsteady  water  and  ice  velocities,  depths,  and  thicknesses. 


REVIEW  OF  ICE  JAM  MODELING 

It  is  convenient  to  review  prior  ice  jam  modeling  in  the  context  of  the  ways  jams 
develop  and  are  classified.  This  section  presents  an  ice  jam  classification  system, 
reviews  past  analyses  of  stationary  jams,  and  briefly  describes  existing  numerical 
models  used  for  predicting  jam  thickness. 

Review  of  ice  jam  classification 

The  International  Association  for  Hydraulic  Research  (IAHR)  Working  Group 
on  River  Ice  Hydraulics  (IAHR  1986)  published  a  state-of-the-art  report  classifying 
the  different  types  of  ice  jams  and  reviewing  techniques  for  their  analysis.  The 
report  defines  ice  jams  as  stationary  accumulations  of  fragmented  or  frazil  ice  that 
restrict  flow.  This  broad  classification  could  include  any  form  of  ice  cover  or  accu¬ 
mulation,  except  for  a  thermally  grown  sheet  ice  cover.  The  classification  system 
distinguishes  ice  jams  by  their  season  of  formation,  dominant  formation  process, 
spatial  extent,  and  state  of  evolution.  It  is  clear  from  the  classification  that  jam  for¬ 
mation,  whichever  type  of  jam  forms,  is  intrinsically  unsteady.  Jams  develop  and 
adjust  in  thickness  and  extent  in  accordance  with  flow  conditions,  ice  availability, 
and  weather  conditions.  Also  clear  from  the  report,  however,  is  that  existing  for¬ 
mulations  of  jams  assume  steady  conditions,  although  unsteady  water  and  ice  move- 
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Figure  5.  Cross  section  of  a  freezeup  jam. 


merit  play  important  parts  in  the  development  of  nearly  all  the  jam  types  in  the 
IAHR  classification. 

Season  of  occurrence 

Ice  jams  are  typically  called  freezeup  or  breakup  jams,  in  accordance  with  the 
season  in  which  they  form.  This  classification  suitably  represents  the  hydrological 
and  meteorological  conditions  prevailing  during  jam  formation. 

Freezeup  jams  form  during  periods  of  subfreezing  air  temperatures  when  frazil 
ice  production  is  great.  Their  composition  is  mainly  frazil  and  broken  pieces  of 
shore  ice,  as  depicted  in  Figure  5.  Frazil  ice  forms  in  areas  of  high  water  velocity 
and  turbulence,  where  heat  loss  from  the  water  surface  is  greatly  increased.  Some 
areas  may  remain  open  and  produce  frazil  throughout  the  winter.  Subfreezing  air 
temperatures  reduce  basin  runoff,  resulting  in  fairly  steady  water  discharge  from 
base  flow.  The  ice  discharge  varies  as  frazil  production  increases  or  freezeup  jams 
form  and  cut  off  the  supply  of  frazil  to  downstream  reaches. 

When  frazil  travels  downstream,  it  agglomerates  as  slush,  which  rises  to  the 
water  surface  and  forms  ice  pans.  The  pans  may  break  upon  passing  through  very 


Figure  6.  Surface  jam  resulting  from  juxtaposition  of  frazil  pans. 
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Figure  7.  Freezeup  jam  following  shoving  and  thickening. 


turbulent  areas  of  flow  or  strengthen  by  freeze-thickening.  Pans  may  slow  at  chan¬ 
nel  constrictions  or  stop  at  downstream  ice  covers.  Stopped  and  juxtaposed  frazil 
pans  form  surface  jams  as  shown  in  Figure  6.  With  continued  transport  of  frazil 
into  a  reach,  the  length  of  a  cover  of  juxtaposed  ice  pans  may  increase  to  the  point 
where  downstream  forces  exceed  the  cover's  strength,  causing  shoving,  collapse, 
and  thickening  of  the  cover.  Figure  7  shows  a  freezeup  jam  formed  by  shoving  and 
thickening  of  ice.  Freezeup  jams  may  strengthen  by  surface  freezing,  which  usually 
causes  them  to  be  thinner  than  breakup  jams  formed  at  similar  flow  rates. 

While  water  discharge  may  be  steady  during  the  formation  and  evolution  of  a 
freezeup  jam,  flow  depth,  ice  velocity,  and  jam  thickness  are  not.  With  continued 
frazil  production,  a  freezeup  jam  may  progress  upstream  with  time,  raising  water 
levels  as  it  progresses.  Figure  8,  for  example,  shows  the  temporal  variation  of  the 
location  of  the  upstream  edge  of  a  freezeup  jam  on  the  Salmon  River  in  Idaho  for 
the  winter  season  1990-91.  The  jam  initiated  on  Julian  Day  61  (30  November)  and 
reached  a  maximum  length  of  approximately  34  km  (21  river  miles).  This  freezeup 
jam  is  an  annual  occurrence  and  consistently  initiates  in  a  deep  pool  at  river-mile 
233.  The  upstream  progression  of  the  freezeup  jam  varies  from  year  to  year  and 
depends  on  air  temperature.  Plotted  in  Figure  8  is  the  daily  mean  air  temperature; 
frazil  ice  is  generated  when  air  temperature  plunges  below  0°C.  Water  depths 
experience  unsteady  variations  because  of  the  shoving  and  thickening  of  the 
freezeup  jam.  Indeed,  field  measurements  show  that  as  the  jam  progresses  through 
a  reach,  the  water  level  increases  by  approximately  2  m. 

Breakup  jams  occur  during  periods  of  relatively  warm  weather  and  are  typified 
by  periods  of  increased  runoff.  The  runoff  results  from  snowmelt,  rain,  or  ground- 
water  release.  For  these  jams,  water  discharge  usually  is  highly  unsteady  with  surges 
being  possible  owing  to  the  failure  and  reformation  of  jams  as  breakup  progresses 
downstream.  Many  anecdotal  accounts  of  the  highly  unsteady  nature  of  breakup 
jams  have  been  presented  (Moberley  and  Cameron  1929,  p.  151).  Moberley  was  the 
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Figure  8.  Temporal  variation  of 
air  temperature  and  location  of 
upstream  edge  ofafreezeupjam 
on  the  Salmon  River ,  Idaho , 
1990-91. 


Factor  of  the  Hudson's  Bay  Co.  post  at  Fort  McMurray  in  Alberta  and  described  the 
following  events  on  the  morning  of  20  April  1875: 

...The  winter  of  1874-75  was  a  bitter  one,  with  deep  snow  and  never  a  thaw  until 
April.  On  the  2nd  or  3rd  of  that  month,  however,  a  further  heavy  fall  of  snow  was 
followed  by  a  sudden  rise  in  temperature.  The  change  of  weather  and  the  weight  of 
melting  snow  caused  the  ice  for  the  85  mile  stretch  of  rapids  above  the  fort  [Fort 
McMurray]  to  breakup,  and  it  came  down  the  Athabasca  with  terrific  force.  On  strik¬ 
ing  the  turn  of  the  stream  at  the  post  it  blocked  the  river  and  drove  the  ice  2  miles  up 
the  Clearwater  [a  major  tributary]  in  piles  40  to  50  feet  high.  In  less  than  an  hour  the 
water  rose  57  feet,  flooding  the  whole  flat  and  mowing  down  trees,  some  3  ft.  in  diam¬ 
eter,  like  grass.... 

As  its  name  implies,  a  breakup  jam  consists  of  pieces  of  broken  sheet  ice,  refro¬ 
zen  frazil  ice,  and  brash  or  slush  ice.  Figure  9  depicts  a  fairly  typical  breakup  jam. 
Freeze-bonding  of  ice  pieces  is  usually  negligible  for  a  breakup  jam,  because  the 
above-freezing  air  in  which  they  form  inhibits  it.  Breakup  jams  typically  form  at 
reaches  where  the  downstream  progression  of  a  run  of  moving  ice  (or  breaking 
front)  slows  because  of  reductions  in  channel  slope,  or  width,  or  where  it  encoun¬ 
ters  resistant  portions  of  ice  cover,  such  as  locations  of  freezeup  accumulations. 
The  severity  of  flooding  during  breakup  jams  depends  on  many  factors,  such  as 
initial  ice  cover  thickness  and  strength,  characteristics  of  the  runoff  hydrograph, 
and,  relatedly,  weather.  Gradually  warming  weather  with  no  rain,  for  instance. 


8 


Figure  10.  Mid-winter  breakup  jam  on  the  Kennebec  River  in  Maine. 


mildly  increases  runoff  discharge  and  decreases  ice  strength,  often  resulting  in  less 
severe  jamming  and  flooding. 

Conversely,  mid-winter  jams  happen  with  the  onset  of  mid-winter  thaws  (usu¬ 
ally  in  January  for  the  northern  U.S.,  but  also  in  early  February,  and  late  Decem¬ 
ber).  Rain  and  snowmelt  runoff  on  impermeable  frozen  ground  can  result  in  very 
steep  increases  in  river  discharge  that  break  up  relatively  strong  ice  covers.  While 
ice  covers  in  early  to  mid-winter  are  typically  not  as  thick  as  they  might  be  in  late 
winter  or  early  spring,  mid-winter  jams  can  produce  severe  flooding.  Also,  since 
the  weather  systems  bringing  mild  mid-winter  weather  are  usually  short-lived, 
and  are  soon  followed  by  frigid  weather,  these  jams  may  remain  in  place,  consoli¬ 
dating  as  virtually  a  monolithic  mass  of  ice.  This  sets  the  stage  for  additional  prob¬ 
lems  later  during  the  normal  breakup  period.  Figure  10  shows  a  mid- win  ter  breakup 
jam  on  the  Kennebec  River  in  Maine.  It  formed  during  January  1996  at  a  peak  dis¬ 
charge  of  approximately  2000  m3/s.  Once  the  flow  receded  to  the  river's  normal 
winter  levels  of  200  to  300  m3/s,  the  jam  grounded  in  many  locations  and  contin¬ 
ued  to  cause  increased  water  levels  upstream. 

As  a  further  example  of  the  highly  unsteady  nature  of  breakup  jams.  Figure  11 
shows  the  stage  hydrograph  for  a  gauging  station  on  the  St.  John  River  in  northern 
Maine,  where  breakup  jams  are  an  annual  occurrence.  Superimposed  on  the  gen¬ 
eral  rise  in  river  stage  are  several  short-duration  peaks  attributable  to  ice  jam  for¬ 
mation  and  failure.  Additional  instrumentation  installed  at  the  gauge  site  identi¬ 
fied  the  initial  time  of  failure  of  the  sheet  ice  cover  as  0610  on  22  April  1994,  which 
corresponds  to  the  stage  drop  following  the  first  large  peak.  Subsequent  peaks  are 
most  likely  ascribable  to  reformation  and  failure  of  additional  jams  downstream 
from  the  gauge  or  from  surges  due  to  failure  of  jams  upstream  from  the  gauge. 
Field  observations  have  shown  that  jams  do  form  at  locations  approximately  1, 4, 8, 
and  15  km  downstream  from  the  gauge. 


Figure  11.  Stage  record  for 
USGS  station  no.  0101000  on 
the  St.  John  River  in  Maine  dur¬ 
ing  1 994  breakup. 


Dominant  formation  process 

Ice  jams  are  also  classified  in  terms  of  the  processes  dominating  their  formation. 
The  important  processes  are  transport  and  deposition  of  ice,  congestion  of  ice,  and 
shoving  and  thickening  of  ice. 

Jams  known  as  "hanging  dams"  are  local,  thick  accumulations  of  frazil  ice.  Frazil 
ice  arriving  at  the  upstream  edge  of  an  existing  ice  cover  may  submerge  and  travel 
beneath  the  cover  until  reaching  an  area  of  lower  water  velocity,  where  it  deposits 
on  the  underside  of  the  cover.  These  formations,  or  hanging  dams,  may  continue  to 
develop  during  the  entire  winter  season,  growing  up  to  thicknesses  of  20  m.  Large 
hanging  dams  may  constrict  flow  and  block  the  downstream  passage  of  ice,  occa¬ 
sionally  initiating  more  severe  breakup  jams. 

A  congestion  or  surface  jam  forms  when  ice  transport  along  the  water  surface  is 
reduced  by  shore  ice  growth,  transverse  floating  objects,  or  channel  constrictions, 
such  as  bridges.  This  initiation  method  is  typical  of  both  freezeup  and  breakup 
jams.  Water  flow  is  fairly  steady  for  congestion  jams,  which,  initially  at  least,  are 
relatively  thin  (i.e.,  single  floe  thickness). 

A  submergence  or  frontal  progression  type  of  jam  typically  forms  at  some  trans¬ 
verse  floating  barrier,  such  as  an  intact  ice  cover.  Arriving  ice  floes  tip  and  sub¬ 
merge,  but  come  to  rest  almost  immediately  at  the  upstream  edge  of  the  obstruc¬ 
tion,  causing  the  jam  to  progress  upstream.  This  type  of  jam,  often  called  a  narrow 
channel  jam,  primarily  forms  during  the  freezeup  of  contacting  frazil  pans. 

Given  a  steady  ice  supply  from  upstream,  both  the  congestion  and  frontal  pro¬ 
gression  type  of  jams  may  grow  to  a  point  beyond  which  they  are  not  able  to  with¬ 
stand  the  increased  stream  wise  load  exerted  by  additional  ice  conveyed  to  the  jam 
or  by  increased  discharge.  Then,  shoving  and  thickening  occur.  A  jam  formed  by 
shoving  and  thickening  is  sometimes  termed  a  wide  channel  jam.  It  will  remain  in 
place  as  long  as  the  downstream  acting  forces  of  water  shear  and  gravity  (and  pos¬ 
sibly  wind  shear)  can  be  resisted  by  the  jam's  strength  and  resistive  shear  stresses 
acting  between  the  ice  and  the  banks.  If  the  load  (water  drag  and  ice  weight) 
becomes  too  great,  the  jam  must  thicken  to  increase  its  resistance  to  downstream 
movement.  Thickening  increases  both  the  jam's  strength  and  the  shear  stress 
between  the  jam  and  the  bank. 

Once  a  jam  collapses,  thickens,  and  possibly  moves  downstream,  jam  reforma¬ 
tion  becomes  highly  unsteady.  The  shear  stress  exerted  on  the  jam  underside 
reduces  as  it  begins  moving,  which  in  turn  increases  water  velocity  and  reduces 
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depth.  Ice  acceleration  produces  an  ice  momentum  that  must  be  overcome  to  arrest 
ice  motion  at  a  location  downstream  where  thickening  takes  place.  During  shoving 
and  thickening  of  accumulating  ice,  water  depth  and  velocity,  ice  thickness,  and  ice 
velocity  are  all  interrelated  and  vary  with  distance  and  time.  It  is  this  unsteady 
nature  of  shoving  and  thickening  that  this  study  addresses. 

Spatial  extent 

The  horizontal  and  vertical  extent  of  a  jam  are  also  used  to  classify  jam  type.  In 
plan,  jams  are  either  partial  or  complete,  according  to  their  extent  across  a  river.  A 
partial  jam  means  that  a  portion  of  the  river  width  remains  as  an  intact  ice  sheet,  or 
one  channel  around  an  island  jams  while  the  other  remains  clear.  In  vertical  sec¬ 
tion,  jams  are  classified  as  floating  or  grounded.  Grounding,  when  ice  extends  to 
the  channel  bed,  takes  place  quite  often  near  the  riverbanks,  at  shallow  areas  such 
as  bars  or  crossings,  and  near  the  toe  region  of  a  thick  jam.  Grounded  jams  usually 
result  from  very  unsteady  water  and  ice  flows.  They  severely  limit  water  flow,  greatly 
increasing  water  levels.  Water  may  flow  as  seepage  through  grounded  accumula¬ 
tions  or  even  over  the  top  of  a  jam.  Little  is  known  about  the  mechanism  of  ground¬ 
ing  or  the  permeability  of  grounded  jams.  Floating  jams  are  more  common  and 
easier  to  analyze,  though  they  may  become  partially  grounded  when  river  flow 
recedes.  Most  analyses  to  date  assume  floating  jams,  whose  flotation  follows 
hydrostatic  pressure  law. 

State  of  evolution 

The  final  classification  category  is  that  describing  the  state  of  jam  evolution: 
steady-state,  evolving  in  time,  or  evolving  up-channel.  An  evolving  jam  continues 
to  be  subject  to  unsteady  flow  rates,  ice  discharges,  or  changes  in  other  ice  variables 
(such  as  strength).  A  breakup  jam,  already  formed  and  undergoing  shoving  and 
thickening,  will  continue  evolving  with  nonuniform  thickness,  depth,  and  water 
velocities.  A  freezeup  jam  may  experience  fairly  steady  flow  rates,  but  frazil  ice 
production  and  transport  may  cause  it  to  shove  and  thicken  with  time.  Figure  12 
shows  ice  jam  evolution  with  time.  As  the  jam  thickens  and  progresses  upstream, 
water  levels  rise  and  velocities  decrease.  Whatever  the  final  water  surface  level  and 
jam  thickness  profile  might  look  like,  the  ice  thickness  and  velocity  are  not  steady 
as  the  jam  develops. 

If  conditions  do  become  steady,  uniform,  and  stable,  a  jam  may  have  an  equilib¬ 
rium  section.  Strictly  speaking,  an  equilibrium  section  is  uniform  only  in  a  reach- 
averaged  sense.  Figure  12d  shows  that  the  thickness  and  depth  are  nearly  constant 
in  the  longitudinal  direction.  Moreover,  the  bed  slope,  water-surface  slope,  and 
energy  slope  are  equivalent  in  the  equilibrium  section.  These  conditions  of  unifor¬ 
mity  were  assumed  by  ice  researchers  when  formulating  the  first  analyses  of  forces 
exerted  on  a  stationary  jam  (e.g.,  Pariset  and  Hauser  1961,  Uzuner  and  Kennedy 
1976). 

Analysis  of  stationary  jams 

A  major  advance  in  addressing  the  effects  of  ice  jams  on  water  levels,  and  in 
estimating  jam  thicknesses,  was  the  realization  that  a  floating  jam  could  be  likened 
to  a  granular  material  contained  between  two  parallel  walls.  The  behavior  of  a 
granular  material  is  influenced  by  the  forces  exerted  upon  it  and  its  material  prop¬ 
erties.  As  the  length  of  a  jam  increases,  these  forces  increase,  as  do  stresses  within 
the  jam. 
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The  first  analysis  was  done  by  the  Canadian  R J.  Kennedy  (1958),  who  was  inter¬ 
ested  in  the  forces  exerted  on  a  boom  by  a  pulpwood  jam.  To  determine  those  forces, 
he  developed  the  following  relation 

B-  +  xiB-2k0XF  =  0  (1) 


where 

B  =  width  of  the  holding  area 

Tj  =  shear  stress  on  the  underside  due  to  water  flow 
fc0  =  coefficient  of  lateral  thrust 

X  -  coefficient  of  friction  of  the  pulpwood  against  the  shore  boundary 
F  =  force  per  unit  width  acting  in  the  downstream  direction  (x-direction). 


If  there  exists  a  forced  exerted  by  the  water  impinging  on  the  upstream  end  of  the 
accumulation,  then  eq  1  can  be  integrated  to  obtain  the  force  at  any  location  x  (mea¬ 
sured  from  the  upstream  end  of  the  accumulation),  i.e. 
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Equation  2  shows  that,  as  the  length  of  the  accumulation  increases,  the  force 
level  reaches  an  asymptotic  value  that  is  related  only  to  the  water  drag  and  cover 
width  (friction  and  lateral  stress  coefficients  being  constant).  Kennedy's  formula¬ 
tion,  however,  neglected  the  streamwise  component  of  the  weight  of  the  pulp  wood 
and  assumed  that  the  depth  (and  therefore  water  shear)  was  constant.  Consequently, 
the  force  through  the  jam  would  be  independent  of  the  thickness  and  bulk  density 
of  the  pulpwood  accumulation. 

Berdennikov  (1964)  investigated  the  forces  exerted  on  an  ice  boom  retaining  an 
ice  accumulation.  While  he  initially  identified  the  weight  component  of  the  ice  mass 
parallel  to  the  water  surface  as  one  of  the  forces  to  be  considered,  he  dismissed  this 
force  and  the  hydrodynamic  pressure  against  the  leading  edge  of  the  cover  /x  as 
being  so  small  as  to  be  negligible.  His  expression  for  the  normal  stress  in  the  ice 
field  gx,  assuming  that  Gx  =  0  at  x  =  0,  is 


Gv  = - — 

2M0ti 


1  -  exp 


-2  Xk0x 
B 


(3) 


Pariset  and  Hausser  (1961),  then  Pariset  et  al.  (1966),  advanced  Kennedy's  for¬ 
mulation  by  including  the  streamwise  components  of  the  weight  of  the  cover  and 
an  assumed  "cohesive"  stress  acting  between  the  ice  and  the  banks.  Summation  of 
the  forces  acting  on  the  ice  cover  (Fig.  13)  gives 


dFB  +  2(xcT|  +  Xk0F)dx  =  (xj  +/3  )Bdx  (4) 

where 

F  =  force  per  unit  width  acting  in  the  downstream  direction 
xc  =  cohesive  stress  per  unit  area  at  the  banks 
r|  =  coverthickness 

/3  =  downstream  component  of  the  weight  of  the  cover  per  unit  area  and  B ,  fc0, 
X,  and  Xj  are  defined  as  above,  so  that 

/3  =  wns  (5) 

where 

Sj  =  specific  gravity  of  ice 
p  =  density  of  water 
g  =  acceleration  due  to  gravity 
S  =  slope. 


Figure  13.  Forces  acting  on  an  ice  cover.  (After  Pariset  et  al.  1966.) 
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Because  their  formulation  is  for  a  uniform  channel,  bed  slope,  water  surface  slope, 
and  energy  slope  are  taken  as  equal.  Integration  of  eq  4  results  in  an  expression  for 
the  longitudinal  force  per  unit  width  as  a  function  of  the  distance  from  the 
upstream  edge  of  the  cover,  similar  to  eq  2,  i.e. 


F  =  - 
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wher efi  is  as  defined  above.  Equation  6  is  based  on  an  assumed  "equilibrium  thick¬ 
ness"  of  an  ice  jam  over  steady,  uniform  flow.  Pariset  et  al.  also  suggest  definitions 
for  narrow  and  wide  jams.  When  the  term  within  the  large  brackets  of  eq  6  (the 
multiplier  of  the  exponential  term)  is  negative,  the  longitudinal  force  F  is  a  maxi¬ 
mum  at  the  upstream  edge  of  the  cover  (x  =  0).  This  is  the  case  for  so-called  narrow 
jams.  As  the  cover  progresses  upstream,  the  downstream  thrust  is  resisted  by  shear 
stress  at  the  banks,  which  grows  faster  than  the  additional  hydrodynamic  forces 
exerted  on  the  jam.  Conversely,  when  the  term  within  the  brackets  is  positive,  the 
longitudinal  force  F  grows  with  distance  downstream  from  the  upstream  edge  of 
the  jam,  asymptotically  approaching  a  maximum  as  the  distance  grows  very  large. 
This  maximum  longitudinal  force  acting  through  a  wide  jam  is 
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Pariset  et  al.  recognized  that  this  maximum  force  (or  sum  of  external  forces)  is 
resisted  by  the  strength  of  the  accumulated  ice,  which  is  assumed  to  behave  as  a 
granular  material.  If  Fmax  exceeds  jam  strength,  the  jam  fails  and  must  thicken  until 
there  is  a  balance  between  the  external  forces  and  jam  strength.  They  likened  the 
maximum  strength  of  the  ice  jam  to  that  of  a  granular  material  under  complete 
mobilization  of  the  passive  pressure  resistance,  i.e. 

KpSiPg{l -Si)^-=tan2j^  + 1  jsjfig(l -S;)^-  (8) 

where  Kp  is  a  passive  pressure  coefficient  and  <j>  is  the  angle  of  internal  resistance  of 
the  accumulated  broken  ice,  and  is  commonly  taken  as  the  angle  of  repose  for  granu¬ 
lar  materials.  Pariset  et  al.  then  equated  jam  internal  stress  to  the  sum  of  the  stresses 
exerted  by  the  external  forces.  In  doing  so,  they  introduced  the  coefficient  (i,  where 


\i  =  k0KpX  (9) 

which  combines  the  ice  properties  into  one  coefficient.  The  stress  balance  results  in 
an  equation  relating  jam  thickness  and  stresses  exerted  against  the  jam: 

•  (10) 

The  shear  stress  Tj  in  eq  10  is  expressible  as 
n2 

(11) 

where  u  is  water  velocity  beneath  the  cover  and  C  is  the  Chezy  coefficient.  Pariset 
et  al.  assumed  that  values  of  C  are  essentially  equivalent  for  the  ice  surface  and  the 
bed. 
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(12) 


The  weight  of  the  cover  in  the  downstream  direction  can  be  defined  as 


/3=SiPgil-2 


C  R] 


H 


where  RH  is  the  total  hydraulic  radius  beneath  the  ice.  Substituting  eq  11  and  12 
into  10,  and  dividing  by  H2  (H  =  the  open  water  depth  just  upstream  of  the  ice 
cover),  renders  eq  10  dimensionless,  i.e. 
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(13) 


Equation  13  can  be  used  to  predict  ice  thickness  for  the  wide-jam  case.  Pariset  et  al. 
recognized  that  freeze-bond  forces  are  of  lesser  importance  for  thicker  jams,  and 
are  not  important  during  breakup  conditions,  because  jam  resistance  is  dominated 
by  gravitational  effects,  i.e.,  ice  weight.  They  developed  a  dimensionless  stability 
parameter  X  to  relate  jam  thickness  to  upstream  open-water  conditions: 


(14) 


where  Q  is  water  discharge  and  wQis  bulk  open  water  velocity  for  uniform  flow  far 
upstream  of  the  jam.  Values  of  X  are  plotted  in  Figure  14  for  Sj  =  0.92  and  |ll  =  1.28, 
which  are  values  adopted  by  Pariset  et  al.  for  jams  in  the  St.  Lawrence  River.  The 
figure  indicates  jams  as  being  stable  (inside  the  bell  curve)  or  unstable  (outside  the 
curve).  The  curve  is  useful,  but  it  does  not  enable  direct  calculation  of  ice  thickness. 
Moreover,  it  assumes  equilibrium  conditions,  i.e.,  steady,  uniform  flow  of  water 
and  uniform  ice  thickness. 

Uzuner  and  Kennedy  (1976)  presented  a  detailed  formulation  of  the  time- 
dependent  differential  equations  describing  the  force  equilibrium  in  a  static,  float¬ 
ing  ice  jam.  Their  formulation  is 
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Figure  14.  Dimensionless  stability 
parameter.  (After  Pariset  et  al.  1966.) 


15 


(15) 


^CTx  ^  ^  (nTxy ) ' “  Ti  C0S(Si  -  siP^n  sin(e  +  a)  =  0 

where 

Gx  =  normal  stress  in  the  streamwise  direction 
TXy  =  shear  stress  at  the  banks 
Tj  =  shear  on  the  underside  of  the  cover 
0  =  slope  of  the  bed 
(0  +  a)=  slope  of  the  water  surface. 

Uzuner  and  Kennedy  expressed  ox  and  Txy  as  functions  of  the  average  vertical 
stress  gz  within  the  cover 

<^=2  siPS(l  -  Si X1  -  V )  cos(9  +  «)  =  Yell  (16) 

where  p  is  jam  porosity  and  ye  is  the  equivalent  unit  weight  of  the  jam.  The  Rankine 
and  Mohr-Coulomb  stress  theories  for  granular  materials  give 

ax  =  Kpoz  (17) 

and 

^xy=^'o(7z+Q  (18) 

where 

Kp  =  passive  pressure  coefficient 
CD  =  shear  stress  coefficient 
Q  =  assumed  cohesive  intercept. 


Substitution  of  eq  16  through  18  into  15,  integration  of  the  modified  equation,  eq 
15,  then  normalization  using  xQ  =  x/hn  and  r|0  =  r| /hn,  yields 
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where  T. 
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2  Kpye 

fc3=  C°/ln 

KpB 

(22) 

n~Us0J  • 

(23) 

Also 

S0  =  bed  slope  (sin  0) 

/b  =  Darcy- Weisbach  resistance  factor  of  the  bed 

qn  =  unit  discharge  at  a  location  upstream  where  ice  does  not  affect  the  flow. 
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The  Uzuner  and  Kennedy  formulation  of  the  force  balance  associated  with  gradu¬ 
ally  varied,  unsteady  water  flow  was  too  complex  for  a  general  solution.  They  did 
show  that,  for  a  condition  of  quasi-steady  jam  formation  in  which  the  jam  progresses 
upstream  at  a  constant  rate,  the  unsteady  water  flow  equations  are  constant  with 
time.  For  the  section  of  the  jam  considered  to  be  in  equilibrium,  eq  19  can  easily  be 
solved  for  r|0  using  the  quadratic  formula  by  setting  the  term  on  the  left-hand  side 
to  zero. 

The  formulations  proposed  by  Pariset  et  al.  and  Uzuner  and  Kennedy  compose 
the  basis  for  most  subsequent  analyses  of  static  ice  jams.  Beltaos  (1983),  most  nota¬ 
bly,  adapted  the  formulations  for  wide-river  jams  and  expressed  flow  depth 
beneath  a  jam  h  as 
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(24) 


where  fQ  is  a  composite  value  of  the  Darcy- Weisbach  resistance  coefficients  for  the 
bed  and  the  jam  underside.  Solution  of  eq  10  for  jam  thickness,  assuming  that  cohe¬ 
sion  is  negligible,  and  that/3  is  given  by  eq  5,  yields 
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where /j  is  a  Darcy- Weisbach  resistance  factor  for  flow  along  the  jam  underside. 
Beltaos  also  presented  field  data,  consisting  of  thickness  measurements  for  several 
ice  jams  that  had  refrozen  in  place.  Using  eq  25,  he  back-calculated  values  of  ji  and 
found  them  to  range  from  0.6  to  3.5,  with  these  upper  and  lower  limits  obtained  for 
conditions  of  considerable  uncertainty.  If  the  two  extreme  values  are  excluded,  his 
data  show  consistently  that  (i  =  0.8  to  1.3.  Beltaos  found,  on  average,  that  \i  =1.2, 
which  is  in  good  agreement  with  the  value  of  1.28  suggested  by  Pariset  et  al. 

Numerical  modeling 

Several  numerical  models  of  jams  have  been  developed.  They  assume  that  a 
balance  exists  between  forces  acting  on  the  jam,  predict  equilibrium  jam  thickness, 
and  estimate  jam  effects  on  water  levels.  Existing  open-water  models  for  steady 
and  unsteady  flow  simulations  have  been  adapted  by  the  use  of  equations  similar 
to  eq  25  to  provide  estimates  of  ice  jam  conditions.  Other  models  have  been  devel¬ 
oped  using  steady  or  unsteady  water  flow  and  equilibrium  (uniform)  or 
nonequilibrium  thickness. 

HEC-2,  the  step-backwater  program  developed  by  the  U.S.  Army  Corps  of  Engi¬ 
neers,  was  modified  to  include  an  ice  cover  (HEC  1979).  In  its  initial  version,  the 
cover,  or  jam,  was  treated  simply  as  a  boundary,  floating  at  hydrostatic  pressure, 
that  provides  an  additional  resistance  to  flow  at  the  water  surface.  The  cover  is 
taken  as  being  static,  with  the  user  of  the  program  inputting  values  of  cover  thick¬ 
ness,  roughness,  location,  specific  gravity,  and  a  value  of  |i.  The  program  calculates 
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the  ice-affected  water  levels  with  the  input  configuration  of  the  cover  and  also 
determines  if  the  cover  is  stable  according  to  the  stability  parameter  proposed  by 
Pariset  et  al.  To  generate  an  equilibrium  thickness  profile,  many  iterations  are  nec¬ 
essary  with  modifications  made  to  the  ice  thickness  and  roughness  values. 

DWOPER,  an  unsteady  flow  forecasting  model  developed  by  the  U.S.  National 
Weather  Service,  was  also  adapted  to  investigate  the  effects  of  ice  covers  on  water 
levels.  Daly  and  Ashton  (1983)  modified  the  St.  Venant  equations  describing  the 
unsteady  water  flow  to  include  the  frictional  resistance  of  the  ice  to  the  water  flow. 
They  concentrated  on  running  steady  water  discharges  with  a  stationary  ice  cover 
and  then  instantaneously  removing  the  cover,  simulating  a  complete  cover  failure 
and  passage  downstream.  The  ensuing  transients  increased  with  increasing  bed 
slopes  as  would  be  expected.  Their  work  did  not  include  nonuniform  covers  or 
jams  and  they  pointed  out  the  necessity  of  including  ice  motion  in  developing  a 
truly  unsteady  ice  jam  model. 

Flato  and  Gerard  (1986)  and  Flato  (1987)  applied  their  model,  ICEJAM,  to  pro¬ 
duce  ice  jam  profiles  for  a  range  of  steady-flow  discharges.  As  it  uses  a  form  of  the 
differential  equation  describing  the  balance  of  forces  on  the  ice  cover,  similar  to 
that  of  eq  19,  it  can  be  used  to  describe  the  complete  thickness  profile  even  if  there 
is  not  an  equilibrium  section.  The  input  data  necessary  to  run  ICEJAM  include 
water  discharge,  ice  jam  characteristics  (bulk  specific  gravity,  angle  of  internal 
resistance,  and  porosity),  channel  data,  roughness  of  the  bed  and  ice,  and  initial 
estimates  of  water  depth  and  jam  thickness.  The  model  first  calculates  the  normal 
depth  (under  ice)  profile  based  on  the  initial  estimates  of  ice  thickness.  It  then  solves 
the  ice  force  balance  equation  in  a  forward-difference  mode,  stepping  downstream 
from  the  upstream  end  of  the  jam.  The  hydraulic  conditions  are  then  modified  for 
these  new  ice  thicknesses  by  means  of  the  standard  step-backwater  calculation  tech¬ 
nique  moving  in  an  upstream  direction.  Iterations  of  the  ice  and  water  calculations 
continue  until  an  acceptable  tolerance  is  met.  Adjustments  are  made  in  the  ice  thick¬ 
ness  at  the  toe  of  the  jam  in  relation  to  a  prescribed  ice  erosion  velocity.  The  model 
produced  reasonable  and  stable  results  when  a  damping  factor  of  1  /$  was  applied 
to  the  calculated  corrections  for  ice  thickness. 

RIVJAM  is  a  model  developed  by  Beltaos  (1993).  It  is  based  on  a  similar  model 
proposed  by  Beltaos  and  Wong  (1986).  Both  models  use  a  steady  water  discharge 
and  include  the  seepage  flow  through  the  jam  in  an  attempt  to  better  define  ice 
thickness  near  the  toe  of  the  jam,  which  may  be  grounded.  RIVJAM  solves  first- 
order  differential  equations  for  the  water  depth  beneath  the  cover  and  the  ice  thick¬ 
ness.  It  does  so  by  means  of  a  predictor-corrector  scheme,  and  the  solution  proce¬ 
dure  may  progress  in  an  upstream  or  downstream  direction.  Beltaos  (1993)  showed 
that  RIVJAM  was  able  to  reproduce  ice  thickness  profiles  for  a  variety  of 
nonequilibrium  and  potentially  grounded  jams  quite  well,  with  appropriate  choices 
of  several  model  parameters.  The  most  tenuous  of  these  appears  to  be  the  seepage 
coefficient,  which  is  similar  in  concept  to  hydraulic  conductivity  (with  units  of 
length /time)  for  high  Reynold's  number  flows.  The  model,  however,  does  not 
include  the  unsteady  movement  of  the  ice  cover  and  thus  cannot  include  the  effects 
of  ice  momentum,  which  would  be  important  in  cases  of  grounded  jams. 

A  utility  program  was  developed  by  Wuebben  et  al.  (1995)  for  use  with  HEC-2  to 
simplify  calculation  of  ice-affected  water  levels.  The  program,  dubbed  ICETHK, 
uses  standard  output  variables  from  a  HEC-2  simulation  and  calculates  ice  thick¬ 
ness  based  on  an  equation  similar  to  eq  25.  There  are  several  calculation  options, 
such  as  width  smoothing,  ice  thickness  smoothing,  and  overbank  ice  and  rough- 
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ness  coefficient  assignment  as  a  function  of  thickness  based  on  the  data  of 
Nezhihovskiy  (1964).  The  program  then  automatically  updates  the  HEC-2  input 
file  to  reflect  these  new  values  of  ice  thickness  and  roughness.  Iterations  continue 
until  a  specified  tolerance  is  met.  Considerable  judgment  is  necessary  in  the  jam- 
toe  area,  where  ice  thickness  conditions  cannot  be  expressed  adequately  by  the 
equilibrium  thickness  as  provided  by  eq  25. 

Lai  and  Shen  (1991)  developed  the  jam  model  RICE,  which  is  intended  to  simu¬ 
late  unsteady  conditions  of  water  flow,  water  temperature,  ice  concentration,  and 
thermal  growth  and  decay  of  ice.  In  their  model,  ice  travels  downstream  at  the 
water  velocity  until  it  reaches  some  location  where  a  jam  forms,  by  either  ice-piece 
juxtaposition  or  the  narrow-jam  or  wide-jam  accumulation  modes.  The  wide-jam 
mode  is  taken  to  be  governed  by  the  ice  force  balance  equation  proposed  by  Pariset 
et  al.  Lai  and  Shen  did  recognize  that  as  progression  (by  shoving)  is  taking  place, 
there  is  a  simultaneous  change  in  the  flow  hydraulics.  They  take  care  of  flow  changes 
by  solving  the  equilibrium  thickness  and  step-backwater  equations  simultaneously 
in  the  reach  where  the  jam  is  thickening.  The  RICE  model  has  been  used  success¬ 
fully  in  simulations  of  ice  conditions  on  the  St.  Lawrence,  Niagara,  Ohio,  and  Yel¬ 
low  rivers,  though,  like  other  models,  it  requires  significant  calibration  to  match 
field  data. 

Tsai  et  al.  (1988)  developed  a  jam  model  to  investigate  ice  transport  in  rivers  and 
ice  jam  initiation.  They  used  a  one-dimensional  numerical  scheme  for  solving  the 
ice  transport  equations,  i.e.,  conservation  of  ice  momentum,  ice  mass,  and  ice  area. 
The  equations  are  solved  in  a  Lagrangian  form,  where  the  trajectories  of  ice  ele¬ 
ments  at  fixed  Eulerian  grid  points  at  the  beginning  of  a  time  step  are  traced  on  the 
x-t  plane.  Values  of  ice  variables  are  then  interpolated  back  to  the  grid  points  at 
the  end  of  the  time  step.  The  de  Saint  Venant  equations  for  unsteady  water  flow  are 
solved  using  a  four-point  implicit  finite-difference  scheme.  The  ice  transport  and 
water  flow  equations  are  loosely  coupled  by  first  solving  the  water  flow  equations 
and  then  the  ice  transport  equations  based  on  the  new  values  of  the  water  flow 
variables. 

Shen  et  al.  (1990)  elaborated  further  aspects  of  this  model,  examining  the  vari¬ 
ous  plausible  constitutive  relationships  possible  for  describing  the  internal  stresses 
and  bank  shear.  For  example,  they  describe  a  rapid  flow  regime  as  one  in  which  the 
ice  concentration  is  low  and  interaction  between  ice  particles  is  minimal.  Commen- 
surately,  they  characterize  a  slow  flow  regime  as  one  in  which  higher  (multi-layer) 
ice  concentrations  typically  form,  and  where  internal  resistance  is  attributable  to 
prolonged  interaction  of  contacting  particles.  Their  expressions  for  the  streamwise 
stress  ox  and  the  stress  normal  to  the  bank  Txy  are  equivalent  to  those  for  passive 
pressure  resistance,  as  described  by  Pariset  et  al.  The  authors  state  the  model 
appears  to  adequately  describe  the  time  and  location  of  jam  initiation  in  river  chan¬ 
nels,  but  that  more  research  is  necessary  to  improve  the  constitutive  laws. 

Summary 

While  considerable  progress  has  been  made  in  modeling  the  unsteady  flow 
associated  with  stationary  ice  jams,  the  unsteady  aspects  of  ice  movement  have 
not  been  adequately  addressed.  Most  models  treat  shoving  and  thickening  as  an 
instantaneous  phenomenon,  with  no  consideration  for  the  effects  of  ice  momen¬ 
tum  on  the  resulting  jam  thickness  and  profile.  Physics  and  field  observations  sug¬ 
gest  that  ice  momentum  should  substantially  affect  jam  thickness. 

The  following  sections  describe  laboratory  and  numerical  experiments  aimed  at 
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evaluating  the  effects  of  ice  momentum  on  jam  thickness  and  profile.  The  numeri¬ 
cal  model  used  for  this  purpose  is  a  significant  advance  on  prior  models  in  so  far 
that  it  includes  ice  momentum  and  directly  couples  ice  and  water  motion. 


LABORATORY  EXPERIMENTS 
Introduction 

As  indicated  in  the  last  section,  the  literature  on  ice  jams  contains  no  studies 
describing  how  jam  shoving  and  thickening  occur,  or  generally  evaluating  the 
importance  of  ice  momentum  in  jam  development.  The  laboratory  experiments 
conducted  here  provide  the  first  diagnostic  information  demonstrating  the  impor¬ 
tance  of  ice  momentum. 

All  prior  studies,  certainly  those  that  do  not  include  ice  motion,  treat  shoving 
and  thickening  as  an  instantaneous  process.  When  the  forces  exerted  in  the  down¬ 
stream  direction  on  a  jam  reach  the  level  of  the  passive  pressure  resistance  of  the 
jam,  prior  formulations  let  the  jam  simply  thicken.  No  mention  is  made  of  the  time 
required  for  thickening  to  take  place,  or  where  the  ice  mass  required  for  the  thick¬ 
ening  originates.  Equilibrium  thickness  theory  (e.g.,  see  Uzuner  and  Kennedy  1976) 
carries  with  it  many  assumptions,  including  steady,  uniform  flow  and  a  stationary 
ice  cover.  Certainly,  when  a  jam  fails,  it  violates  the  latter  assumption,  which  in 
turn  violates  the  steady  and  uniform  flow  condition,  because  ice  movement  influ¬ 
ences  water  flow.  Once  an  ice  jam  comes  into  motion,  the  shear  stress  on  the  under¬ 
side  of  the  jam  is  reduced,  because  it  is  a  function  of  the  difference  between  water 
and  ice  velocities.  Furthermore,  the  principal  assumption  used  for  describing  the 
compressive  stress  state  of  ice  jams  diminishes  in  validity  once  a  jam  fails.  The 
Mohr-Coulomb  theory  has  been  used  with  great  success  in  describing  the  com¬ 
pressive  strength  of  granular  materials,  such  as  ice  rubble  in  a  jam,  under  various 
states  of  stress.  Once  failure  begins,  however,  the  material  undergoes  changes  in 
stress  levels  that  are  not  well  handled  using  this  theory.  As  well  as  thickness  and 
velocity  changes,  other  jam  characteristics,  such  as  porosity  or  even  ice-piece  size 
or  shape,  may  change. 

To  model  shoving  and  thickening,  the  principal  effect  of  ice  momentum,  it  is 
necessary  to  know  how  the  process  occurs.  Though  numerous  ice  jams  and  their 
failures  have  been  observed  for  a  wide  variety  of  situations  in  the  field,  observa¬ 
tions  are  typically  limited  to  the  surface  of  the  cover  from  the  perspective  of  the 
shoreline.  Even  when  jam  failure  and  reformation  are  observed  from  the  air,  prac¬ 
tical  limitations  (i.e.,  altitude  and  sight  distance)  render  the  observations  reach- 
averaged  at  best.  The  highly  unsteady  nature  of  most  breakup  jams  reduces  oppor¬ 
tunities  for  direct  measurements  of  jam  properties.  Only  in  the  rare  incidence  where 
a  jam  formed  and  refroze  in  place,  following  a  reduction  in  water  flow  and  a  return 
of  lower  air  temperatures,  might  this  be  done.  While  these  few  cases  may  provide 
useful  data  on  jam  thickness  profiles,  other  items  of  interest,  such  as  ice  velocity 
and  local  water  discharge  at  the  time  of  jamming,  remain  unknown.  A  final  note 
concerning  jam  observations  is  that,  while  the  date  of  breakup  ice  runs  and  jam¬ 
ming  in  the  northern  U.S.  might  average  10  March  (near  equal  amounts  of  daylight 
and  darkness  daily),  about  80%  of  ice  runs  and  jams  take  place  during  darkness. 

To  qualitatively  examine  the  importance  of  ice  momentum  on  jam  processes,  a 
laboratory  study  was  undertaken  to  simulate  the  shoving  and  thickening  of  a  fail¬ 
ing  ice  cover.  Of  particular  interest  are  the  timing  and  mechanics  of  the  process. 
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The  study  used  a  laboratory  flume  and  a  jam  formed  of  ice  pieces  and  plastic  beads. 
The  jam  was  destabilized  by  means  of  flow  increases.  The  qualitative  observations 
provided  the  insights  necessary  to  numerically  model  unsteady  jam  formation.  A 
brief  series  of  quantitative  experiments  was  also  undertaken  to  address  the  appli¬ 
cability  of  equilibrium  thickness  theory  for  determining  jam  thickness,  following  a 
shoving  and  thickening  event. 

Experimental  setup 

The  experiments  were  undertaken  in  three  different  flumes  using  real  and  plas¬ 
tic  ice.  Two  of  the  flumes  are  refrigerated.  One  flume  has  a  fixed  slope  bed  (S0  = 
0.0033),  is  1.22  m  wide  with  a  working  depth  of  0.61  m,  and  is  22.9  m  long.  The 
second  flume  is  tiltable,  36.6  m  long  with  a  working  section  of  1.22  m  wide  by  0.61 
m  deep  (Fig.  15).  Both  flumes  are  housed  in  refrigerated  rooms  in  which  air  tem¬ 
peratures  can  be  regulated  to  -23°C.  A  selection  of  pumps  delivers  water  from  large 
sumps  to  each  flume,  resulting  in  non-recirculating  flow.  The  sumps  contain  chiller 
coils  to  further  reduce  water  temperature.  The  quantitative  experiments  were  con¬ 
ducted  using  the  tiltable  flume  and  a  third,  unrefrigerated  flume.  This  latter  flume 
has  a  fixed  horizontal  bed,  is  7.3  m  long,  and  has  a  working  section  that  is  0.92  m 
wide  by  0.92  m  deep.  It  has  a  large  sump  and  can  be  operated  in  either  the  fully 
recirculating  or  "once  through"  mode. 

The  visualization  experiments  entailed  initially  forming  ice  covers  made  of  a 
single  layer  of  ice  or  plastic  beads.  Then,  water  discharge  was  increased  to  destabi¬ 
lize  the  cover  and  induce  shoving  and  thickening.  For  the  experiments  with  real 
ice,  ice  pieces  were  formed  from  a  thin  sheet  (about  15  mm)  grown  in  the  flume  at 
very  low  flow.  The  randomly  shaped  pieces  were  approximately  80-120  mm  in 
their  longest  dimension.  Air  temperature  was  increased  to  about  0°C,  and  the  flow 
was  increased  until  the  ice  pieces  collected  as  a  single-layer  accumulation,  or  jam, 
held  in  place  by  a  screen  at  the  downstream  end  of  the  flume.  The  flow  rate  was 
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Figure  16.  Plastic  beads  used  to  simulate  ice. 

then  further  increased  in  steps  until  the  jam  destabilized  and  shoving  and  thicken¬ 
ing  took  place.  Step  increases  in  flow  rate  were  varied  from  10  to  100%  of  the  initial 
flow  rate. 

Tests  were  also  conducted  using  plastic  beads  formed  from  extruded  polyethyl¬ 
ene  strands  that  are  chopped  to  produce  uniform  pieces  approximately  cylindrical 
in  shape  with  a  length  and  diameter  of  3  mm  (Fig.  16).  The  specific  gravity  of  the 
plastic  is  0.925.  The  beads  form  accumulations  with  a  porosity  of  0.40,  very  similar 
to  natural  ice,  and  have  a  dry  angle  of  repose  of  36°.  When  new  beads  are  added  to 
the  water,  they  exhibit  some  surface  tension  effects  (nonwetting),  but  after  a  few 
days  they  become  fully  wetted.  The  beads  were  used  for  the  visualization  experi¬ 
ments  in  the  tiltable  flume,  as  well  as  for  the  quantitative  testing,  because  of  their 
constant  and  uniform  properties,  as  opposed  to  real  ice.  For  the  visualization  tests, 
a  uniform  layer  of  beads  was  spread  over  the  water  surface  at  a  very  low  flow.  The 
flow  rate  was  then  increased  to  the  starting  flow  level  for  the  test  and  the  bead 
cover  was  allowed  to  consolidate,  typically  resulting  in  thicknesses  of  one  to  two 
beads.  Then,  similar  to  the  real  ice  tests,  the  flow  rate  was  increased  in  steps  until 
shoving  and  thickening  took  place. 

The  quantitative  experiments  were  conducted  in  two  series  of  tests.  The  first 
were  conducted  in  the  unrefrigerated  horizontal-bed  flume  to  determine  impor¬ 
tant  physical  parameters  characterizing  accumulations  of  beads.  In  these  experi¬ 
ments,  a  stable  bead  jam  was  formed  at  a  variety  of  flow  rates,  with  detailed  mea¬ 
surements  taken  of  water  surface  slopes,  depths,  accumulation  thicknesses,  and 
velocity  profiles.  From  these  measurements,  the  values  of  the  composite,  bed,  and 
jam  underside  friction  factors  were  calculated,  as  was  the  value  of  (i  for  the  beads. 
The  second  series  of  tests  was  conducted  using  the  tiltable-bed  flume  in  a  warm 
environment.  Similar  to  the  visualization  experiments,  a  bead  jam  was  formed  and 
disturbed  by  flow-rate  increases.  Measurements  were  made  of  jam  thickness, 
extent,  water  velocity,  water  surface  slope,  and  depth.  These  measurements, 
together  with  the  values  determined  for  |i  and  friction  factors,  allowed  the  jam 
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thickness,  following  shoving  and  thickening,  to  be  compared  with  that  predicted 
by  equilibrium  thickness  theory. 

Observations  of  shoving  and  thickening 

The  jams  failed  in  two  general  ways  with  both  the  real  ice  and  plastic  beads 
simulating  ice.  The  type  of  failure  was  related  to  the  initial  flow  rate,  jam  thickness, 
and  the  relative  increase  in  flow  rate  over  the  initial  value. 

The  first  mode  of  jam  failure  is  called  here  "progressive  jam  failure."  It  was  the 
dominant  mode  for  low  initial  water  discharge,  relative  to  the  discharge  needed  to 
fail  the  entire  cover,  and  when  the  relative  increase  in  discharge  was  less  than  50%. 
As  the  discharge  was  increased,  the  water  level  rose  first  at  the  upstream  end  of  the 
flume  as  the  discharge  wave  traveled  its  length.  The  consequent  increase  in  shear 
stress  on  the  underside  of  the  cover  caused  minor  consolidation  at  the  upstream 
end  of  the  accumulation  and  some  undertuming  and  transport  of  individual  pieces. 

Small  ridges  of  local  thickening  formed  near  the  jam's  upstream  end.  The  ridges 
developed  rather  slowly,  with  the  jam  upstream  moving  into  the  ridge,  while  the 
portion  of  the  jam  downstream  remained  motionless.  The  increased  thickness  and 
roughness  of  the  ridge  further  increased  shear  stress  on  the  jam's  underside  in  the 
vicinity  of  the  ridge,  initiating  additional  small  failures  and  ridge-building  events 
further  downstream.  When  new  ridges  formed  downstream,  the  activity  at  the 
upstream  ridges  slowed  or  ceased  and  the  entire  jam  above  the  most  downstream 
ridge  began  moving. 

With  time,  this  progression  of  ridge-forming,  herein  termed  the  "shoving  front," 
advanced  to  the  jam's  downstream  end  at  the  screen.  At  that  point,  the  entire  jam 
was  in  motion.  The  jam  thickness  subsequently  increased  at  the  screen  until  the 
thickening  jam's  strength  could  resist  the  downstream-acting  forces.  Although  ridge¬ 
building  resulted  in  minor  local  thickening  as  it  progressed  downstream,  the  major 
thickening  occurred  at  the  screen,  then  progressed  back  upstream.  As  thickening 
progressed  upstream,  the  portion  of  the  jam  downstream  of  the  "thickening  front" 
became  stationary. 

Meanwhile,  the  jam  upstream  was  still  moving.  Figure  17  provides  an  idealized 
picture  of  the  movements  of  the  shoving  front  and  the  thickening  front  during  the 
progressive  failure  of  a  jam.  Ice  velocities  during  jam  failure  and  thickening  were 
very  much  less  than  the  bulk  water  velocity. 
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Figure  17.  Progressive  jam  failure. 
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The  second  mode  of  jam  failure  is  called  here  "complete  jam  failure/'  It  was  the 
dominant  mode  for  relatively  high  initial  discharge,  for  which  the  entire  jam  was 
close  to  a  condition  of  instability.  Discharge  increases  simply  overwhelmed  the 
entire  jam.  As  discharge  increased,  water  level  rose,  and  a  wave  of  water  traveled 
the  length  of  the  flume.  In  these  events,  however,  the  entire  ice  jam  (which  was  up 
to  30  m  long)  mobilized  en  masse  and  failed  at  the  downstream  screen.  The  thick¬ 
ening  front  then  progressed  upstream  from  the  screen.  During  some  tests  with  very 
large  discharge  increases,  small  ridges  formed  elsewhere  within  the  jam,  but  the 
major  thickening  took  place  as  the  thickening  front  swept  back  upstream.  The  ice 
velocities  for  this  type  of  failure  were  noticeably  higher  than  for  the  progressive 
failure  type,  yet  still  only  ranged  up  to  25%  of  the  bulk  water  velocity. 

A  few  tests  were  also  conducted  in  which  the  discharge  was  increased  and  held 
constant  for  a  short  period.  These  tests  were  conducted  for  initial  discharges  and 
discharge  increases  previously  identified  as  causing  a  progressive  jam  failure.  In 
these  tests,  the  shoving  front  was  allowed  to  progress  about  halfway  down  the 
flume,  then  the  discharge  was  reduced  to  its  original  rate.  As  the  discharge 
receded,  the  entire  jam  (which  was  moving  upstream  of  the  shoving  front)  stopped 
en  masse.  This  left  an  ice  accumulation  that  was  nearly  a  single  layer  thick  in  the 
downstream  reaches,  but  slightly  thicker  upstream.  The  discharge  was  then 
increased  to  the  higher  value.  In  all  cases,  the  jam  upstream  of  where  the  shoving 
front  had  previously  progressed  mobilized  en  masse.  The  shoving  front  then  con¬ 
tinued  its  progression  downstream  as  if  the  drop  and  subsequent  increase  in  dis¬ 
charge  had  never  happened. 

While  the  failure  modes  observed  for  the  ice  pieces  and  plastic  beads  were  gen¬ 
erally  similar,  there  were  a  few  differences.  As  expected,  the  real  ice  was  more 
angular  and  thus  had  a  higher  angle  of  internal  resistance.  Consequently,  the  accu¬ 
mulations  of  real  ice  were  more  resistant  to  increases  in  downstream  load.  A  com¬ 
plication  for  the  tests  using  ice  was  water-temperature  regulation.  Depending  on 
air  temperature,  an  accumulation  may  melt  or  it  may  further  increase  its  strength 
owing  to  freeze-bonding  of  contacting  pieces.  The  ice  pieces  were  also  much  larger 
than  the  plastic  beads,  potentially  violating  assumptions  that  their  behavior  could 
be  treated  using  continuum  or  particulate  theory. 

An  interesting  finding  from  the  tests  was  that  small  increases  in  discharge  do 
not  necessarily  result  in  shoving  and  thickening.  Sometimes,  two  or  three  small 
steps  in  discharge  were  required  to  destabilize  a  jam.  This  was  especially  true  for 
the  tests  using  ice,  which  involved  a  single  layer  of  ice  pieces  with  a  rather  high 
piece  aspect  ratio  (L/q0),  whereas  the  bead  experiments  involved  a  layer  thickness 
of  between  one  and  two  bead  diameters.  The  high  aspect  ratio  for  the  ice  pieces 
invalidates  assumptions  of  continuum  theory  for  treating  jam  strength  behavior. 
The  initial  discharge  used  for  a  test  likely  was  substantially  less  than  the  discharge 
needed  to  destabilize  the  jam.  For  each  bead  experiment,  once  a  shoving  and  thick¬ 
ening  event  was  completed,  the  water  discharge  was  further  increased  to  start  a 
second  shoving  and  thickening  event.  Usually,  multiple  discharge  steps  were 
required  for  this,  especially  when  the  initial  event  had  led  to  complete  jam  failure. 
Evidently,  the  collapsed  jam  had  thickened  to  an  extent  much  greater  than  the  equi¬ 
librium  thickness  estimated  from  the  existing  formulations  (e.g.,  Uzuner  and 
Kennedy  1976,  Beltaos  1983).  This  preliminary  finding  strongly  suggests  the 
importance  of  ice  momentum  in  determining  jam  thickness.  It  is  this  finding  that 
prompted  further  experiments  to  compare  thicknesses  after  shoving  to  those  cal¬ 
culated  using  equilibrium  theory. 
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Equilibrium  thickness  evaluations 

To  compare  thicknesses  after  failure  and  thicknesses  predicted  from  equilibrium 
theory,  the  strength  and  hydraulic  roughness  properties  of  bead  accumulations  had 
to  be  determined.  A  useful  formulation  of  static  equilibrium  thickness  is 
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in  which  u  is  average  velocity  and/j  is  the  Darcy- Weisbach  friction  factor  for  the 
flow  along  the  ice  cover  underside.  The  overall  strength  coefficient  for  the  ice  p  is  a 
combination  of  several  material  properties,  i.e. 


H  =  k0XKp(l-p) 


(27) 


where 

k$  =  lateral  stress  coefficient  (the  fraction  of  the  longitudinal  acting  force  that 
is  directed  normal  to  the  banks) 

X  =  friction  coefficient  for  ice  sliding  against  ice  at  a  shear  boundary 
p  =  accumulation  porosity 
iCp  =  Rankine  passive  pressure  coefficient 


with  0  being  the  angle  of  internal  resistance,  which  is  commonly  assumed  to  be 
equal  to  the  dry  angle  of  repose  of  a  granular  material. 

Calculation  of  equilibrium  thickness  using  eq  26  requires  knowledge  about  the 
average  velocity  and  energy  slope  of  the  flow,  as  well  as  about  the  specific  gravity 
of  the  ice,  accumulation  porosity,  friction  factor  for  flow  along  the  accumulation, 
and  the  overall  strength  coefficient  p.  For  jam  equilibrium,  bed  slope,  water  surface 
slope,  and  energy  slope  are  taken  as  being  equal.  While  values  of  angle  of  internal 
resistance  <j)  and  porosity  p  had  been  directly  measured  for  the  beads,  values  of 
lateral  stress  coefficient  k0  or  the  friction  coefficient  X  are  unknown.  Therefore,  p 
cannot  be  readily  calculated  from  eq  27.  Instead,  values  were  back-calculated  using 
eq  26,  but  this  approach  involves/j  as  an  additional  unknown.  However,  the  Darcy- 
Weisbach  definition  of  friction  slope  for  the  ice-affected  layer  of  flow 
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combined  with  eq  26  leads  to  an  alternate  form  of  eq  26 
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For  this  equation,  only  values  of  slope  and  hydraulic  radius  of  the  ice-affected  flow 
area  are  needed  to  evaluate  r)eq. 

Multi-layer  bead  accumulations,  or  jams,  were  allowed  to  form  in  the  flume  at 
different  levels  of  steady  discharge.  Detailed  slope  measurements  and  velocity  pro¬ 
files  were  obtained  in  the  region  where  the  accumulation  thickness  appeared  uni¬ 
form.  Figure  18  presents  an  example  of  a  measured  velocity  profile  and  the  fitted 
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Figure  18.  Measured  velocity 
profile  with  fitted  log-law  equa¬ 
tions  for  the  ice  and  bed-affected 
areas. 


log-law  profiles  calculated  for  the  bed-  and  ice-affected  portions  of  the  flow  area. 
By  assuming  the  line  of  zero  shear  stress  to  be  equivalent  to  the  intersection  of  the 
two  computed  profiles,  the  depth  of  the  ice-affected  flow  area  could  be  determined 
and  the  hydraulic  radius  R j  of  that  region  calculated.  The  value  of  |i  was  then  calcu¬ 
lated  by  back  substitution  into  eq  30,  using  the  measured  value  of  accumulation 
thickness  in  the  equilibrium  reach.  The  value  of  ]±  was  found  to  average  0.75  for  the 
three  steady  flow  discharges  considered.  Though  this  value  is  at  the  low  end  of  the 
range  of  0.8-1. 3  reported  by  Beltaos  (1995),  it  reflects  the  effect  of  the  difference 
between  the  shape  of  the  beads  and  natural  ice  rubble.  The  beads  are  uniform  in 
size  and  approximately  cylindrical  in  shape,  thereby  having  an  angle  of  repose  less 
than  that  of  natural  ice.  The  uniform  shape  of  the  beads  causes  their  accumulations 
to  deform  more  easily  under  stress  and  results  in  a  lower  and  thus  1 1  value  than 
is  the  case  for  natural  ice  rubble  in  a  jam. 

With  a  known  average  value  of  jli  determined  for  the  beads,  further  experiments 
could  proceed  using  the  tiltable-bed  flume  without  refrigeration.  Fifteen  experi¬ 
ments  were  conducted  under  a  variety  of  initial  discharges  and  discharge  increases, 
as  reported  by  Zufelt  (1992).  Each  experiment  was  run  in  a  manner  similar  to  the 
visualization  series,  for  which  a  bead  cover  was  allowed  to  form  at  a  low  discharge 
and  the  flow  then  increased  in  steps  until  a  shoving  and  thickening  event  took 
place.  The  energy  slope  was  assumed  to  be  equivalent  to  the  water  surface  slope, 
which  was  calculated  from  measurements  of  water  surface  elevation  along  the  flume. 
Slope  was  plotted  against  discharge,  and  a  linear  relation  was  obtained  for  the 
conditions  before  and  after  failure.  Though  a  power  relationship  is  to  be  expected 
between  slope  and  discharge,  the  linear  fit  was  adequate  for  the  limited  range  of 
discharges  used.  The  two  relations  also  confirm  an  increase  in  jam  roughness  fol¬ 
lowing  failure  and  thickening,  as  evidenced  by  the  higher  slopes  following  failure 
(Fig.  19). 

The  variations  in  slope  for  similar  discharges  before  failure  in  Figure  19  reflect 
slight  variations  in  the  configuration  of  the  bead  jam  (thickness,  extent,  etc.) 
between  successive  tests.  The  slope-discharge  relation  was  used  for  further  calcu¬ 
lations.  The  composite  Darcy- Weisbach  friction  factor  fQ  was  calculated  for  each 
experiment  as 
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Figure  19.  Slope  vs.  discharge,  showing  in-  Figure  20.  Calculated  equilibrium  jam 
crease  in  slope  after  jam  failure.  thickness  vs.  discharge. 


in  which  u  is  average  flow  velocity  and  R  is  the  composite  or  total  hydraulic  radius. 
An  average  value  of  fQ  was  calculated  and  used  for  the  conditions  before  and  after 
failure.  The  detailed  velocity  profiles  from  each  test  provided  data  on  the  ratio  of 
the  bed-affected  to  ice-affected  hydraulic  radii.  Then,  using  the  Sabaneev  equa¬ 
tions 


fo  = 


fi  +/b 
2 


(32) 


and 
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it  is  possible  to  determine  ratios  of  f  to and/j  to  fQ.  Average  values  of  these  ratios 
for  the  conditions  before  and  after  failure  were  calculated  and  used  to  calculate  the 
equilibrium  cover  thickness  expected  for  each  experiment.  Figure  20  is  a  plot  of  the 
equilibrium  thickness  vs.  discharge  for  this  series  of  experiments,  using  average 
values  of  fQ  andf  //G. 

Figure  21  follows  the  changes  in  thickness  experienced  with  step  changes  in 
discharge  for  one  of  the  experiments.  The  initial  thickness  of  a  bead  jam  is  slightly 
greater  than  that  predicted  using  the  equilibrium  theory  (eq  30)  and  was  between 
one  and  two  beads  thick.  Two  step  increases  in  discharge  were  necessary  before  a 
shoving  and  thickening  event  occurred,  with  the  resulting  thickness  again  slightly 
greater  than  that  expected  by  theory.  This  first  failure  was  a  progressive  jam  failure 
with  very  small  ice  velocities.  Two  more  step  increases  in  discharge  were  necessary 
to  again  cause  jam  failure.  The  second  failure  was  a  complete  jam  failure,  with  the 
whole  cover  mobilizing  en  masse  and  thickening  taking  place  initially  at  the  down¬ 
stream  screen.  The  final  thickness  was  significantly  greater  than  the  equilibrium 
value  for  that  flow  level. 

The  jam  failures  in  each  experiment  were  identified  as  either  progressive  or  com¬ 
plete  jam  failures.  The  final  thickness  following  failure  was  plotted  in  Figure  22 
against  the  equilibrium  thickness  predicted  using  eq  30.  The  progressive  jam  fail¬ 
ures  resulted  in  accumulation  thicknesses  that  plotted  on  or  very  slightly  above  the 
equality  line  in  the  figure.  The  complete  jam  failures,  however,  exceeded  the  equi¬ 
librium  thickness  in  every  case,  often  by  a  significant  amount. 
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Figure  21.  Change  in  jam  thickness  with 
step  increases  in  discharge.  A  progressive 
and  a  complete  cover  failure  are  shown. 


t,  Predicted  (mm) 


Figure  22.  Measured  vs.  predicted  jam 
thickness  following  progressive  and  com¬ 
plete  jam  failures. 


Discussion 

The  laboratory  experiments  were  conducted  for  two  reasons:  to  obtain  prelimi¬ 
nary  insights  into  the  shoving  and  thickening  process,  and  to  examine  the  applica¬ 
bility  of  equilibrium  thickness  theory  for  determining  jam  thickness  in  unsteady 
water  and  ice  flow  situations.  The  degree  of  ice  movement  and  its  effects  on  the 
water  flow  appeared  to  depend  not  only  on  the  material  properties  of  the  particu¬ 
late  material  composing  the  jam,  but  also  on  the  initial  discharge  and  subsequent 
discharge  increase.  Two  failure  modes  were  identified:  progressive  and  complete 
jam  failure.  Progressive  jam  failure  took  place  at  lower  initial  discharge  and  lower 
discharge  increases  relative  to  the  discharge  needed  to  completely  destabilize  the 
jam.  It  can  be  characterized  as  the  smooth  movement  of  a  shoving  front  that  travels 
downstream  through  the  jam  to  the  downstream  end,  causing  minor  consolidation 
and  thickening.  On  reaching  the  downstream  screen,  a  thickening  front  moves  back 
upstream,  resulting  in  a  new,  greater  jam  thickness  able  to  withstand  the  higher 
discharge  rate.  For  progressive  jam  failures,  the  final  jam  thickness  is  very  close  to, 
yet  slightly  greater  than,  the  thickness  predicted  using  equilibrium  jam  theory  as 
represented  in  eq  30. 

Complete  jam  failure,  on  the  other  hand,  occurred  for  initial  discharges  close  to 
the  discharge  necessary  to  completely  destabilize  the  cover.  It  can  be  characterized 
by  the  absence  of  (or  instantaneous  passage  of)  the  shoving  front.  The  entire  jam  is 
moved  downstream  en  masse,  failing  and  thickening  at  the  downstream  screen. 
The  final  thickness  of  jam  is  greater,  sometimes  significantly  so,  than  the  equilib¬ 
rium  thickness  predicted  using  eq  30. 

The  experiments,  particularly  those  producing  complete  jam  failure,  point  to 
the  importance  of  ice  momentum  in  determining  the  thickness  resulting  from  the 
arrest  of  a  moving  ice  jam.  The  experiments  also  indicate  that,  especially  for  the 
progressive  jam  failures,  the  time  for  a  shoving  and  thickening  event  to  occur  can 
be  quite  significant  since  the  ice  velocities  are  quite  low.  Shoving  and  thickening 
are  not  instantaneous  and  should  probably  not  be  treated  as  such.  While  the  pro¬ 
gressive  jam  failure  mode  may  not  result  in  significantly  greater  ice  thickness  than 
that  predicted  using  equilibrium  jam  theory,  there  arise  conditions  of  nonuniformity 
of  ice  jam  thickness,  velocity,  and  depth  that  could  significantly  alter  the  very  forces 
that  determine  water  and  ice  flow. 

Additionally,  the  experiments  reveal  the  important  interactions  of  ice  movement 
and  water  flow  on  one  another.  They  show,  therefore,  the  necessity  of  using  a  fully 
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coupled  numerical  model  of  ice  jam  formation.  Such  a  model  includes  both  ice 
velocity  and  the  effects  of  ice  momentum  on  the  force  balance,  utilizing  the  full 
conservation  of  mass  and  momentum  equations  for  the  ice. 


FORMULATION 
Introduction  and  assumptions 

Formulated  here  are  the  one-dimensional,  unsteady  flow  equations  for  water 
and  ice.  The  equations,  derived  in  integral  form,  are  based  on  the  conservation  of 
mass  and  momentum  for  water  and  ice  flow  during  jam  formation  and  breakup. 
The  integral  equations  are  then  discretized  as  finite-difference  equations,  approxi¬ 
mating  the  conservation  laws  in  their  integral  form.  The  equations  are  expressed  in 
terms  of  four  dependent  variables  that  fully  describe  the  flow,  as  shown  in  Figure 
23,  namely  the  velocities  of  the  ice  cover  and  under-ice  water  flow  (u  and  u,  respec¬ 
tively),  the  ice  cover  thickness  (r|),  and  the  under-ice  water  depth  (d).  All  four  vari¬ 
ables  are  functions  of  time  and  space.  The  equations  are  derived  first  in  a  general 
form,  and  then  simplified  in  accordance  with  the  assumptions  listed  below.  Addi¬ 
tional  simplifications  for  certain  flow  conditions  are  discussed  subsequently. 

The  assumptions  made  in  developing  the  equations  for  water  and  ice  flow 
include  the  usual  St.  Venant  assumptions  for  one-dimensional  flows  (e.g.,  Cunge  et 
al.  1980)  are  as  follows: 

•  Flow  is  one-dimensional,  with  velocity  uniform  across  each  cross  section,  and 
water  level  horizontal  at  each  cross  section. 

•  Streamline  curvature  is  small  and  vertical  accelerations  are  negligible,  so  that 
pressure  distribution  is  hydrostatic. 

•  The  effects  of  turbulence  and  boundary  friction  can  be  accounted  for  through 
resistance  laws  identical  to  those  used  for  steady-state  flow. 

•  Average  channel-bed  slope  is  small  so  that  the  cosine  of  the  angle  it  makes 
with  the  horizontal  may  be  taken  as  unity.  (This  assumption  is  valid  for  bed 
slopes  to  about  0.01.) 


d(x,t) 


Figure  23.  Ice  and  ivaterfloiv. 


29 


Other  assumptions  include: 

•  Cross  sections  are  uniform  and  prismatic  in  shape,  ensuring  that  streamline 
curvature  remains  small.  For  this  analysis,  a  rectangular  channel  shape  is  as¬ 
sumed. 

•  All  flow  is  subcritical.  Though  there  may  be  considerable  changes  in  depth 
and  ice  thickness  between  two  cross  sections  (notably  at  the  locations  of  shov¬ 
ing  or  thickening  fronts),  it  is  assumed  that  through  methods  (i.e.,  with  no 
explicit  representation  of  fronts)  suitably  describe  these  fronts. 

•  Ice-piece  properties  remain  constant  (i.e.,  no  heat  transfer,  phase  change,  or 
freeze-bonding  between  ice  pieces). 

•  Jams  are  particulate  continua,  such  that  forces  and  stresses  are  describable 
using  Mohr-Coulomb  stress  theory  and  an  average  value  across  the  cross  sec¬ 
tion. 

•  Jams  float  with  a  constant  bulk  specific  gravity,  do  not  ground  on  the  channel 
bed,  and  are  not  subject  to  significant  motion  or  accelerations  in  the  vertical 
direction. 


Development  of  equations 

The  integral  form  of  the  equations  for  water  and  ice  flow  are  developed  using  a 
control  volume  approach.  The  Cartesian  coordinate  system  used  is  depicted  in  Fig¬ 
ure  24,  in  which  x  denotes  horizontal  distance  along  the  longitudinal  river  axis,  y 
denotes  vertical  distance,  and  z  denotes  transverse  distance  normal  to  the  longitu¬ 
dinal  axis. 


Conservation  of  water  mass 

Conservation  of  water  mass  requires  the  net  inflow  of  water  entering  a  control 
volume  (bounded  by  Xp  x2,  the  bed,  and  the  bottom  of  the  jam  in  Fig.  24)  during  a 
given  period  be  equal  to  the  change  in  water  storage  within  the  control  volume  for 
the  period,  that  is 

Hr.  ,  .  i  X2\ 


1[(P“A)x2  -(ptM)Xi]df+  J  [(pA)ti  ~(pA)h 
h  x\ 


dx  =  0 


(34) 


For  practical  purposes,  water  is  incompressible,  such  that  p  is  constant  and  eq  34 
reduces  to 


;K  _(wA)x, 


dt  + 


?[(A)t2 -(^)ti]rfx  =  °. 


(35) 


Further  simplifications  are  made  subsequently,  such  as  expressing  area  A  in  terms 
of  flow  depth  d(x,t). 


Conservation  of  ice  mass 

The  net  inflow  of  ice  and  pore  water  (between  the  ice  pieces)  into  the  control 
volume,  bounded  by  Xy  x2,  and  the  bottom  and  top  of  the  jam  in  Figure  24,  is  the 
time  integral  of  the  difference  between  the  mass  flow  rates  entering  the  control 
volume  at  x1  and  leaving  the  control  volume  at  x2,  i.e. 
hr 

j[(pi^i[i-p])Xi  +(pi>Asjp)xj  -(piUA[i-p])x2  —  (puM p)X2 


30 


Figure  24.  Longitudinal  and  cross-sectional  views  of  ice  and  water  flow  areas ,  showing 
coordinate  system  used  in  equation  development. 

where 

pj  =  ice  density 
D  =  ice  velocity 

Aj  =  cross-sectional  area  of  the  jam 
p  =  porosity  of  the  jam 
Sj  =  specific  gravity  of  ice. 

The  first  and  third  terms  in  eq  36  represent  the  mass  flux  of  ice,  while  the  second 
and  fourth  terms  represent  the  pore  water.  Pore  water  is  only  contained  in  that 
portion  of  the  ice  area  below  the  phreatic  surface  (SjAj).  The  experiments  of  White 
(1991)  show  that  the  velocity  of  flow  through  a  stationary  frazil  cover  is  negligibly 
small  (10“5  m/s),  resulting  in  negligible  mass  exchange  between  the  pore  water 
and  the  underlying  water  flow.  Hence,  there  is  no  term  for  seepage  flow  through 
the  jam  provided.  Pore  water  is  assumed  to  move  at  the  same  velocity  as  the  ice 
and  since  p*  =  SjP,  the  ice  and  pore-water  terms  may  be  combined.  Setting  eq  36 
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equal  to  the  change  in  ice  storage  over  the  time  interval  and  considering  p  as  con¬ 
stant  results  in 

hr, 

(  a  \  /  A  \  If#.  r  \  f  A  \  f  A  \  II 

(37) 


1  Mi«i)X2  -(uAjSi)^  dt+  j  [(M)t2  -(M)tJ^  =  0  . 

*1  L 


Further  simplifications  are  made  subsequently,  such  as  expressing  A-x  in  terms  of 
jam  thickness  r|(;c,f) . 


Conservation  of  water  momentum 

The  analysis  examines  the  control  volume  for  water  flow  beneath  an  ice  jam 
whose  channel  cross  section  is  prismatic.  Figure  25  depicts  the  forces  acting  on  the 
control  volume  that  is  bounded  by  Xp  x2,  and  the  bed  and  the  bottom  of  the  jam. 

Conservation  of  momentum  in  the  x-direction  requires  that  the  change  of 
momentum  within  the  control  volume  between  times  tx  and  t2  equal  the  sum  of  the 
net  flux  of  momentum  into  the  control  volume  and  the  integral  of  the  external 
forces  acting  on  the  control  volume  during  the  same  period.  The  momentum  inside 
the  control  volume  at  any  instant  is 
*2 

j{pAu)dx  (38) 

*1 

so  the  net  increase  in  momentum  AM  between  times  h  and  t2  is 


AM=  |  (pAw)^  -(pAu)ti 


dx 


(39) 


The  net  momentum  flux  Mf  into  the  control  volume  between  times  and  t2  is 


Mf=  j 

‘i 


(pA«2)  -(pAw2) 


2  J 


dt 


(40) 


The  external  forces  acting  on  the  water  control  volume  include:  hydrostatic  pres¬ 
sure;  gravity  forces  due  to  the  weight  of  the  water,  ice,  and  pore  water;  and  shear 
stress  at  the  bed,  banks,  and  jam  underside.  The  hydrostatic  pressure  forces  acting 
at  sections  Xi  and  x2  are  and  as  depicted  in  Figure  25.  With  the  level  of  the  phreatic 
surface  above  the  bed  denoted  as  D(x),  the  vertical  distance  above  the  bed  as  8(x), 
and  the  local  width  as  fr(8),  for  any  section  x 


Figure  25.  Forces  acting  on  the 
water  control  volume. 
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(41) 


,  d(x) 

Fpi  =  g  \p[D(x)-S\b(x,b)dd  . 

0 

For  a  rectangular  channel,  b  is  constant  in  the  vertical  and  equal  to  the  top  width  B . 
Therefore 


/  d(x) 

Fpi  =g  Jp[D(x)-8]M5  . 

o 

Substitution  of  D(x)  =  d(x)  +  SjTi(jc)  and  integration  gives 
Fpl' =pgB[d2/2  + dsn]  ■ 

Thus,  the  time  integral  for  the  net  pressure  is 

‘kPldt  =  -  Fp'')dt  =  /|[(p/1)X]  -  (p/a)X2  ]dt 


where 


h 


=  B\d2/l+dsiT\*\  . 


(42) 


(43) 


(44) 


(45) 


Two  gravity  forces  act  vertically  on  the  water  control  volume.  The  first  acts  on 
the  bottom  surface  of  the  control  volume  (the  bed).  It  is  attributable  to  the  combined 
weight  of  water,  ice,  and  pore  water  above.  The  second  acts  on  the  upper  surface  of 
the  control  volume  (the  bottom  of  the  jam).  It  is  attributable  to  the  weight  of  the  ice 
and  pore  water  above.  The  horizontal  component  of  the  first  gravity  force  is 

x2  X2 

Fgl  =  !  (p gA  +  pigAi(l-p)  +  pgA^Sodx  =  j pg(A  +  A^ )S0dx  (46) 


*i 


where  SQ  is  bed  slope 


c 

°  dx 


For  the  period  to  t2 

t2  *2*2 

j Fgldt  =  J  \pg(A  +  Aisi)S0dxdt 
f,  Ml 


(47) 


(48) 


The  gravity  force  attributable  to  the  weight  of  ice  and  pore  water  acting  on  the 
upper  surface  of  the  control  volume  (in  the  ^-direction)  is 


Fg2  =  J  [ptfA (1  -  P)s>b  +  P^ASiPSib }dx=  j  [pgAsiSib ] dx 

x,  x1 

where  Sib  is  the  slope  of  the  jam  underside 

dd 


3(yb  +  d) 
Sib= — 


tyb  |  dd 
dx  dx 


_  Q 

°~lh 


(49) 


(50) 


Substituting  eq  50  into  49  and  integrating  over  the  period  /j  to  f2  gives 


h  hx2 

lfg2 .*=  I  1  PgAisi 
fi  Mi 


isi(so-^jdxdt 


(51) 
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The  remaining  forces  to  be  determined  are  the  boundary  resistance  or  shear  forces. 
The  total  shear  stress  produced  by  the  water  flow  (averaged  for  the  control-volume 
reach)  is 


t  =  pgRSf 


(52) 


where  R  is  the  hydraulic  radius  of  the  section  and  Sf  is  the  friction  or  energy  slope 
associated  with  the  water  flow.  The  Darcy- Weisbach  definition  of  friction  slope  is 


sf  = 


_/o»2 

SRg 


(53) 


where  f0  is  the  composite  (bed  and  ice  cover)  Darcy- Weisbach  resistance  factor. 
Equation  53  substituted  into  eq  52  gives 


z_PgRf0u  =  pfQU 

8Rg  8 


(54) 


This  shear  stress  is  a  total  value,  generated  by  the  differences  in  the  velocities  of 
water  flow  relative  to  the  velocities  of  the  other  boundaries  of  the  control  volume. 
It  can  be  split  into  two  parts:  Tb,  the  shear  stress  at  the  bed  and  bank  boundary,  and 
ij,  the  shear  stress  at  the  ice  boundary  Prior  formulations  (e.g.,  Beltaos  1983)  have 
shown  that  the  simple  case  of  a  static  ice  cover  can  be  analyzed  approximately 
using  a  "two-layer  approach,"  separating  the  total  flow  area  into  one  layer  domi¬ 
nated  by  shear  stress  on  the  bed  and  banks,  and  another  layer  dominated  by  shear 


b) 

Figure  26.  Two-layer  approach  designation  of  shear  stress  due  to  water  flow. 
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stress  on  the  ice  cover.  Figure  26  identifies  the  two  layers.  Note  that  a  refers  to  the 
area  of  each  respective  layer,  P  is  the  wetted  perimeter,  and  the  subscripts  b  and  i 
designate  the  bed-  and  ice-affected  variables.  The  dashed  line  indicates  the  line  of 
nominal  zero  shear  or  the  boundary  between  the  two  layers.  The  total  shear  force 
per  unit  length  of  flow  area  is 

tP  =  +/tiPi  •  (55) 

If  the  "two-layer  approach"  is  valid  for  any  value  of  ice  velocity  u,  such  as  depicted 
in  Figure  27,  then  the  shear  stress  at  each  boundary  is  expressible  as 

(56) 

b  8 


and 


T.  =  p/i(»--u)|»-^|  (57) 

1  8 

The  absolute  value  sign  captures  directional  shear.  It  can  be  dropped  provided  that 
stress  direction  is  preserved  in  the  momentum  equation. 

The  "two-layer  approach"  assumes  that  each  layer  can  be  adequately  described 
using  the  Darcy-Weisbach  relationship  for  flow  resistance  and  thus  can  be  related 
to  the  friction  slope  of  the  water  flow,  i.e. 


s  /bM2  fj  (»--»)2 
f  8  Rbg  8Rig 


(58) 
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As  R=  AfP,  eq  58  can  be  rewritten  as 


(59) 

l«bJ  UJU  A  «  J  (59) 

where  and  otj  are  the  bed  and  ice-affected  portions  of  the  flow  area.  The  total 
shear  stress  is  xb  +  X;  and  the  total  shear  force  on  a  unit  length  of  channel  is  xbPb  + 
XjPj ,  which  leads  to  the  following  expression  for  Sf 

s  T  xP  =%Eb±^jPj 

f  P  gR  P gA  p gA  '  (60) 

The  time  integral  of  the  friction  force  is 


I F(dt  =  J  J  PgASfdxdt  =  J  |  [xbPb  +  XiPijrfxdf.  (61) 

fl  ^1  *1*1 

Equations  55  to  57,  rearranged  in  terms  of  the  unknown  variables  u,  v,  d,  and  t|, 
become 


^bPb+tiPj 


P/b*'2pb  .  Pfiju-vfPj 


P/b“2pb  1  ,  fi{u-v)2Pj  _  P/b»2pb  L  Oj_ 

8  [  /b«2pt  _  8  [  «b. 

An  expression  for  Sf  results  when  eq  59  and  62  are  combined 


c  —  +  _  A?**  2  | 


[  /b  fu  «  ;  j  ^ 

The  effects  of  flow  and  ice  velocities  on  TjPj  are  evident  in  Figure  28.  As  ice  veloc¬ 
ity  increases  beyond  water  velocity,  the  shear  stress  caused  by  the  jam's  presence 
reverses  to  the  downstream  direction  and  the  portion  of  the  flow  area  affected  by 
this  force  increases. 

The  full  momentum  equation  for  the  water  flow  can  be  written  as 


£2  £2  £2  £2 

AM  -  Mf  =  j  Fpldt  +  j  Fgldt  -  J  Fg2dt  -  f  Ffdt . 

£]  tl  £j  £a 

Equations  39, 40,  44, 48,  51,  and  61  combined  with  eq  64  give 


Figure  28.  Shear  force  on  the  ice 
jam  underside  vs.  ice  velocity 
D  (rjPj  vs.  v). 
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(65) 


f  [(p«A)t2 -(pwA)t  ]dx- /  (Pi/2a)x  ~(pw2^)x 
X\  h*-  1  2J 

*/[(pIl)Xl -{ph)X2]dt  +  gl  (piA  +  AiS^SJxdt 

*1  fl*l 


1 1  X\ 


hxi  f  hx2 

-g]  \  PM  S0 -^.Uxdt-gj  J pASfdxdt. 


dx  ) 


hxi 


If  p  is  constant,  eq  65  simplifies  to 


?[m)(,  -K ,]*- 1  (“2/1)>i  -("2/1)Xs  *=«f[('i)x, 


?[ 

V 

t ix2  7)s1  ^2x2 

gj  l  AiSi—dxdt  +  gl  l  A(S0-S()dxdt. 
Mi  6X  Mi 


(66) 


Conservation  of  ice  momentum 

The  ice  control  volume  under  analysis  is  bounded  by  Xp  x2,  and  the  bottom  and 
top  of  the  jam.  Figure  29  depicts  the  control  volume  and  the  forces  acting  on  the 
jam.  The  force  associated  with  wind  drag  at  the  air/ ice  interface  is  neglected  in  this 
formulation,  though  it  could  readily  be  included.  The  momentum  of  the  ice  and 
pore  water  in  the  control  volume  at  any  instant  is 

/  (pi«A [1  ~P}+ PvAiSip)  dx.  (67) 

*1 


From  Pi  =  Sjp,  and  regrouping  the  p  and  (1  -  p)  terms,  the  net  increase  in  momentum 
within  the  control  volume  between  t j  and  t2  becomes 


AM  = 


(68) 


The  net  momentum  flux  into  the  control  volume  is 


[(pyA[i-p])+(p^2Asip)]x  -  [(p;u2 A  [l  -  p]) + (pi)2 Assp)]x^  •  (69) 

Once  again,  from  pj  =  SjP,  and  regrouping  the  p  and  (1  -  p)  terms,  the  net  momen¬ 
tum  flux  into  the  control  volume  between  t ^  and  t2  becomes 


h\ 

Mf=i 

<i  L 


(pu2Asi)  -(pv2Asi)x 


dt. 


(70) 


Figure  29.  Forces  acting  on  the 
ice  control  volume. 
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The  external  forces  acting  on  the  ice  and  pore  water  control  volume  include: 
hydrostatic  pressure,  a  gravity  force  due  to  the  ice  and  pore  water  in  x-direction, 
shear  stress  on  the  underside  of  the  jam  from  the  water  flow,  and  shear  stress  at  the 
banks.  The  strength  of  the  jam  also  resists  ice  motion,  provided  stress  acting  through 
the  jam  does  not  exceed  the  maximum  allowable  longitudinal  stress.  That  limiting 
stress  places  the  jam  in  the  Rankine  passive  state  of  stress. 

The  hydrostatic  pressure  forces  acting  at  sections  x^  and  x2  are  fP2  and  Fp2 
respectively.  With  the  level  of  the  phreatic  surface  above  the  bottom  of  the  jam 
denoted  as  SjT|(x),  the  vertical  distance  above  the  bottom  of  the  jam  denoted  as  8(x), 
and  the  local  width  as  b(S),  then  for  any  section  x 

,  Si  T\(X) 

fp2  =g  1 P  [siT|M  ~  5]  b(x,  8)d8.  (71) 

0 


For  a  rectangular  channel,  b  is  constant  in  the  vertical  and  equal  to  the  top  width  B. 
Thus 


,  Si  i\(x) 

Fp2  =g  }p[siTl(x)-8]Bd5. 
0 


(72) 


Integration  of  eq  72  yields 
'  (sm)2 

fP2  =P^L7L- 


(73) 


Consequently,  the  time  integral  of  the  net  pressure  Fp2  is 


-Fp2")dt  =  gf[(pl2)xi-(pl2)xz]dt 


where 


h 


=  B 


Ml 

2 


(74) 


(75) 


A  gravity  force  acts  on  the  bottom  surface  of  the  ice  control  volume  (the  bottom 
of  the  jam)  from  the  weight  of  the  ice  and  pore  water  above.  The  horizontal  compo¬ 
nent  of  this  gravity  force  equals  and  counterbalances  Fg2 


*2, 


*2r 


%=  J  [PigM1-p)Sib+PgAisiPSib]dx=  J  [P$ASiSib M 

z,  x, 


For  the  period  t1  to  t2 


(76) 


h  hx2  (  dd'\ 

l  Fg3*  =  J  f  P8Aisi  so  -  y  fadt.  (77) 

h  h  x\  ^  ' 

The  shear  stress  on  the  underside  of  the  jam  attributable  to  water  flow  is  equal 
and  opposite  to  the  shear  force  of  the  cover  on  the  water,  i.e. 

h 

=  lFhdt-  (78) 

h 

In  the  case  of  the  jam,  however,  only  the  shear  force  along  the  jam  underside,  or 
TjPj,  is  included.  From  eq  62  above 
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(79) 


T  p  p/j(a--u)2Pj 

11  8 

Substituting  for  Ffj  gives 


*2  f2*2  ^2  pfi(w-'V))2Pi 

JFfjrff  =  J  JtiPidxdf  =J  j  ^ — -dxdt. 

h  hxi  hx\ 


8 


(80) 


From  the  assumption  that  ice  jams  can  be  treated  as  a  particulate  continuum,  the 
vertical  and  horizontal  stresses  within  jams  can  be  related  using  Rankine  stress 
theory.  Previous  researchers,  as  discussed  in  the  Review  of  Ice  Jam  Modeling  section, 
give  the  shear  stress  at  the  banks  Txy  as  a  function  of  the  stress  along  the  jam's 
longitudinal  axis  ox.  For  a  granular  material,  the  normal  stress  can  be  expressed 
directly  as  a  function  of  the  vertical  stress  ov,  which  results  from  the  weight  of  ice 
acting  downward  to  the  phreatic  surface  through  the  jam,  and  from  the  buoyancy 
force  acting  upward  in  the  submerged  portion.  Vertical  stress  ascribable  to  ice  weight 
acting  downward  throughout  the  entire  cover  thickness  is 


avi=PigV{1~V)  (81) 

where  y  is  measured  from  the  top  of  the  ice  surface.  The  stress  of  the  water  acting 
downward  in  the  wetted  portion  is 

avw=Ps(y-|l-Si]n)p  (82) 

where  y  is  measured  from  the  top  surface  of  the  ice  and  y  =  (1  -  sj)r\  represents  the 
phreatic  surface.  The  buoyancy  force  acting  upward  in  the  wetted  portion  of  the 
thickness  is 


m' = PS'ly — E1 — Sijn)-  (83) 

From  eq  81  to  83,  the  total  or  effective  vertical  stress  av  is  zero  at  the  upper  and 
lower  surfaces  of  the  jam,  and  ov  reaches  a  maximum  at  the  phreatic  surface.  Over 
the  entire  thickness,  the  average  vertical  stress  is 

^  =  SiP^(1-y)(l-si)y-  (84) 

The  horizontal  or  longitudinal  stress  in  a  granular  material  relates  to  the  average 
vertical  stress  as 


ox=^i^v  (85) 

where  Aq  is  a  coefficient  of  proportionality  reflecting  the  jam's  state  of  stress,  i.e., 
passive,  neutral,  or  active.  For  passive-pressure  loading  conditions,  Aq  typically 
assumes  the  value  of  the  passive  pressure  coefficient  Kp,  which  represents  the  maxi¬ 
mum  or  failure  stress  of  the  material 


where  <J>  is  the  angle  of  internal  resistance  measured  in  degrees.  Equations  81,  82, 
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83,  and  85  combine  to  give  the  force  due  to  normal  pressure  for  a  general  channel 
shape,  i.e. 


F  = 
rPP 


V 


nM 


1  Pi(l-p)hW-8]b(x,8)rf8' 


SiTl(x) 

-  I  P(l- 


p)[siTi(x)  -  5] &(  a:,  5) 


(87) 


For  a  rectangular  channel,  b  is  constant  in  the  vertical  and  equal  to  the  top  width  B. 
Therefore 

FPP' =KpgB{l-p) 

Integration  yields 

Fpp'=KpgB(l-P)sip(l-Si)i  (89) 

The  time  integral  of  this  net  normal  force  in  the  x-di  recti  on  is 


nfar)  ,  SiTl(x) 

1  Sip[ii(*)-8]d8-  j  p[Siri(x)-8]d8 
o  o 


(88) 


kppdt  =  j(Fpp' ■ ~Fpp")dt  =  ^si(l-si)l[(pf3)xi  -(PI3)XJ* 


where 


h=KpB±-(l-p). 


(90) 


(91) 


It  should  be  recognized  that,  by  setting  the  value  of  /q  =  Kp  in  eq  85,  the  jam  is 
considered  to  be  at  its  strength  limit.  A  jam  at  rest  could  probably  experience  Jq 
values  ranging  from  Ka,  the  active  pressure  coefficient,  to  the  passive  pressure 
coefficient.  For  a  cover  in  motion,  Jq  most  likely  depends  on  ice  velocity.  More 
research  is  needed  to  determine  the  proper  values  of  /q  to  be  used  in  different  states 
of  ice  motion. 

An  important  further  force  is  a  shear  stress  TXy,  produced  by  ice  grinding  against 
the  banks  or  channel  sides.  It  relates  directly  to  the  normal  stress  in  the  x-direction 
at  failure,  i.e. 


'Zxy=<5xk0\  =  cvk0\Kp  (92) 

where  k0  is  the  coefficient  of  lateral  thrust  (percentage  of  normal  stress  acting  in  the 
x-direction  that  would  act  in  the  z-  or  cross-channel  direction)  and  X  is  the  sliding 
coefficient  of  ice  on  ice  at  the  banks.  Bank  shear  stress  acts  equally  at  either  bank 
(for  one-dimensional  formulations,  at  least)  and  may  vary  between  section  Xj  and 
x2.  The  force  from  the  shear  at  the  banks  can  be  obtained  from  eq  84  and  92,  i.e. 


fis  =  2  j  k0XKpsiPg(l  -  p)(l  -  Si )  ^-  Ax. 
X\ 

The  time  integral  of  this  force  is 


(93) 
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(94) 


^2  ^2X2 

f  Fisdt  =  gSi  (1  -  S;  )(l  -  p)  J  I  pl4dxdt 

f,  tjX, 

where 

Ii=k0XK?vl2. 


(95) 


Therefore,  the  full  conservation  of  momentum  equation  for  the  ice  cover 
becomes 


AM  —  Mf  —  J  F^2dt  +  J  Fg3dt  +  J  F^dt  +  J  Fppdt  —  J  F[Sdt 
h  h  h  h  h 


(96) 


or 


J  fep^iL  -(sjP^Ai)  1  dx- I  (pu2AjSi)  -(p\)2A«i 

X,L  2  1J  t,L'  /X1  ' 


=  gj[{ph)Xl  -{Ph)X2]dt  +  gj  J  PM 


hxi  f  p. 

+  Hp  8 
fl*l  8 


(u-v)  dxdt  +  gsill-s^j^pls)^  -(p/3)x 

fl 


rff 


-^i(l-si)(l 


>2*2 

■P)S  jphdxdt. 

hxi 


(97) 


With  p  being  a  constant 


/  —  (sjixAi)  ]dx-j  (t)2Asi)  ~h2Asi)  dt 

XlL  1  1J  fjL  /X1  '  /X2_ 

=  gj[{h)Xj  -{h)X2]dt  +  gj  !  Aisifso-^jdxdt  +  j  j  ^-(u-x>)2dxdt 
+^i(l-si)jT(l3)  -(/3)  'bf-2Si(i_si)(i-p)J  jl4dxdt. 

fi  J  Mi 


(98) 


The  integral  relations  given  by  eq  35,  37,  66,  and  98  are  valid  for  a  channel  of 
constant  rectangular  cross  section.  Substituting  rectangular  channel  cross-section 
relationships  for  Sf,  lv  I2,  /3,  and  /4  (i.e.,  A  =  Bd  and  Aj  =  Bq),  then  canceling  out 
common  terms,  produces  the  conservation  of  mass  and  momentum  equations  in 
their  integral  form. 

Conservation  of  water  mass  is 


x?  r 


*+J[W  t2-Wtl 

x\ 


dx  =  0. 


Conservation  of  ice  mass  is 


(99) 


![H)x2  ~H)xIjdf+  /[fa)t2  -Hj 


dx=  0. 


Conservation  of  water  momentum  is 


(100) 
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i  [(«d),2  -Mt,  ]dx+j  (w2rf)X2  -  (”2rf)Xi  dt 


h 

+8 1 

h 


y+<M 


'*2 


U2  „  ) 

—  +  dSir) 

v  A 


*2*2  0^ 

dt-gj  }  rjSj  —  dxdt 
»i*i 


'2*2  ,  ,  Vf2/b«  (B  +  2d) 

-?J  I  dS0dxdt  +  $  \  ^ — — - 

<i*i  <1*1 


1  + 


u  —  v 


A  b 

fb  {B  +  2d){  u 


dxdt  =  0. 


(101) 


Conservation  of  ice  momentum  is 


J  [H)t,  - HL  ldx  + 1  M  -  dt + gSi  j 

X]  1  2  1J  f,L'  ,X2  /Xlj  f, 


V2/ 


^  n2"1 


x2  ^  'X\ 


<2*2  (  M\  ,  J2 

~g\  \  ?{s0-—jdxdt  +  g(l-si)\ 


<1*1 


Kp{l-p)r\2 


'x9 


Kp(1-P)ti 


2^ 


/X, 


if 


* 


1  *2*2  9  j?(l-S;)f2*2  ,  x  7 

- J  J  fi(u-v)2dxdt  +  — - —  J  |  /c0AfCp(l-p)ri2^f  =  0 

8si  »,*!  B  <1  *1 


(102) 


Discretization  of  the  system  of  equations 

The  equations  are  transformed  from  their  integral  form  into  a  differential  form 
to  facilitate  the  discretization  needed  to  proceed  with  numerically  simulating  jams. 
By  assuming  the  dependent  variables  to  be  continuous,  differentiable  functions 
enables  their  expansion  as  Taylor  series.  The  terms  in  the  conservation  of  water 
mass  equation  become  thereby 


/JX  3(d)  A1  d2{d)At2 
(d).  —(d),  h — — — A f  H - ^ - h... 

V  >l2  v  dt  dt2  2 

,  ,,  ,  ,,  d{ud)  d2(ud)  Ax2 

(ud)  =(ud)  h — r — -Ax  +  — - 

V  ;X2  ^  ^  dx  d  2  2 


(103) 


+  .... 


Then,  by  retaining  only  the  first-order  terms  and  assuming  that  the  increments  Ax 
and  At  approach  zero 


lim 

*2 - 


*2 


m. 


Xl 


dtdx 


/[(“Ox,  -Mx, 


*2*2 


3(wd) 


if  =  j  j  -A- — f-ixif. 
*1*1 


dx 


(104) 


Consequently,  the  conservation  of  water  mass  equation  reduces  to 


*2*2 

JJ 

*i*iL 


d(d)  d(ud)‘ 
dt  dx 


dxdt  =  0. 


(105) 


If  this  equation  is  to  hold  everywhere  in  the  (x,t)  plane,  then  it  will  hold  for  an 
infinitely  small  volume,  such  that  the  eq  105  can  be  rewritten  as 
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(106) 


M+^H= o. 

dt  dx 

In  a  similar  fashion,  the  conservation  of  ice  mass  equation  becomes 

M+iM=0. 

dt  dx 


(107) 


The  conservation  of  water  momentum  and  conservation  of  ice  momentum  equa¬ 
tions  contain  a  combination  of  single  and  double  integrals.  By  Taylor  series  expan¬ 
sion  as  above,  the  conservation  of  water  momentum  becomes 


hx2 

1J 

Mi 


- — -  + 
dt  dx 


)(u2d)  ^ 

fd2  ,  ) 

——  +  rfS;T| 

- - -+  sr  — 

\L  J 

dx 


dxdt 


hx2 

-s\\ 

Mi 


\„ddldS  fbAB  +  2d) 

U ,  fi  B 

'u-v~ 

2Y 

n'dx  dS°  8  gB 

{  fb(B  +  2d ) 

u 

J_ 

dxdt  =  0 


(108) 


and  the  conservation  of  ice  momentum  becomes 


hx2 

1/ 

Mi 


a(UT1)  a(A) 


(„2\ 


dt  +  dx 


+8SiTx 


\  2  J 


hx2 

+gf  s 

Mi 


+s(1-Si)|^Mi-p)l 

k0XK  (1  - p)r I2 -T)S0+r]^--£-(u-v)2 


r  2  "i 

\ 

2 

L  J- 

/ 

dxdt 


dx  8^1 


dxdt  =  0. 


(109) 


Under  the  assumption  that  these  equations  also  hold  everywhere  in  the  (x,t) 
plane,  the  conservation  of  water  momentum  equation  is 


W,a(“2l  3 


(  a1 


dt 


dx  +Sdx 


-  +  ds  jT| 


dd  AC 

-gs^n^-gdSo 


8  B 


i +Il  5  -fUziH 

fb(B  +  2d){  u 


~  0. 


(110) 


In  a  similar  manner,  the  conservation  of  ice  momentum  equation  becomes 
9(uti)  ,  3('°2r>)  ,  „  d  f  r|2  ^ 


dt 


dx  -  +  SSid^ 


v2. 


I  gChfil 


\  y 

k0XKp  (l -  p)t]2  -  gr|S0  +  gr\  ^  -  -k-  (u  -  v)2  =  0. 


(Ill) 


The  partial  derivatives  with  combined  dependent  variables  in  the  foregoing  equa¬ 
tions  are  then  expanded  to  separate  the  variables.  The  conservation  of  water  mass 
equation  becomes 


dd  dd  ,  dll 
-  +  u—  +  d—  =  0. 
dt  dx  dx 


(112) 
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dx>  dri  dd  u  vT,  u  xOri 
■^+Ssi^8^+S^Si)K p(1-p)t- 


id-fi) 


k0XKp  (1  -  p)  ti  ~gS0--r~-{u- v)2  =  0. 


Figure  30  shows  the  computational  grid  used  for  the  numerical  simulations  in 
this  study.  The  discretization  of  the  governing  equations  assumes  that  the  values  of 
the  dependent  variables  and  functions  of  the  variables  between  the  computational 
grid  points  can  be  expressed  in  terms  of  the  values  at  the  grid  points.  A  function 
f(x,i)  within  the  space  interval ;  and  (j  +  1)  is  replaced  by  the  following  weighted 
average 


(116a) 

(116b) 

t 


Figure  30.  Computational  grid  used  for  the  numerical  simulations. 
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where  x^<x<  *j+1.  Similarly,  for  the  time  interval  tn^t<tn+1 

f(xvt)=y^ (ii7a) 

and 

/(*j+i <■ f)  =  e-/j+i 1  +  (a  -  e)/j+i  (117b) 

where  0  <  0  <  1  and  o  <  *¥  <  1  are  weighting  factors.  Subscripts  denote  the  x- 
location  while  superscripts  denote  time  level.  By  use  of  these  relations,  the  finite 
difference  equations  that  approximate  the  conservation  laws  expressed  by  the  dif¬ 
ferential  eq  112,  113,  114,  and  115  can  be  developed.  When  =  0.5,  the  resulting 
equations  represent  the  Preissmann  four-point  scheme  of  finite  differences.  In  the 
discretization,  the  derivatives  of  a  function  f(x,t)  are  estimated  as 


dt 


2At 


and 


dx  Ax 

and  the  values  of  the  variables,  or  the  functions  themselves,  are  estimated  as 


(118a) 

(118b) 


f{x,t)  = 


e(/f+1  +^l1)  +  (l-9)(yin+^l) 


(119) 


The  above  rules  of  discretization  can  be  used  to  recast  the  conservation  of  water 
mass  equation  as 


( dftf+dT1)- 

_ 

e^f+i1 

-df+1)+(i-e)(df+a-df)' 

I  u 

2At 

Ax 

+d 


e^f+i1  -  «”+1 ) + (i  -  e)(t<f+1  -u?) 


Ax 


=  0 


where 


_  _  e« +»r1)+(i-e)(^1+<) 


and 


d  = 


e(^  +  ) + (l  -  e)(^f+i + ) 


Similarly,  the  conservation  of  ice  mass  equation  becomes 


K+i1+<+1)-K+i+<) 


2At 


ej  (t!^1  -  rif +1 )  +  (1  -  0j  -  nf ) 

Ax 


+T[ 


f  n+1  n+ 1) 

(uj+i  J 

+(1-90 

i - 

K.-, 

P 

1 

P 

Ax 

=  0 


(120) 


(121) 


(122) 


(123) 
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where  0j  is  a  weighting  factor  for  the  ice  variables 


-  9ihn++i1+'of+1)+(i-ei)(^+i +•»") 


v  = 


(124) 


and 

_  _  ei  (<+Y  +  T|f +1 )  +  (1  -  6i  )(r|”+1  +  r|" ) 

2 

The  conservation  of  water  momentum  equation  transforms  to 

[uf+i +  !(ji+i) _  (Mi+i + uf )  _r 

2At 


-  +  u 


!  /b»  !  /b»  t  /i  (»-•») 

8rf  4B  8d 


e(«f+Y  -«f+1) 

+(i-e)( 

!uf+i  ~uf) 

I9 

«-rff+r 

Ax 

)+(i-e) 

Kw) 

Ax 

+^i 


e 

i  n+l  rc+ll 

l^j+l  _  Tlj  j 

+  (1-0) 

Ax 

-ss0=o. 


(125) 


(126) 


where  u,d,a nd  +  v  are  defined  as  above.  Note  that  the  last  term  in  eq  114  has  been 
expanded  for  clarity  and  that  the  discretization  of  these  terms  containing  squared 
variables  are  weighted  averages  of  the  squared  terms.  They  are  not  squares  of  the 
weighted  averages  of  the  variables.  This  minor  detail  preserves  the  physical  mean¬ 
ing  of  the  terms.  Those  terms  are 


u2=- 


(127) 


and 


M’C-MT 


+(i-e) 


(128) 


The  conservation  of  ice  momentum  equation  transforms  to 


(uffi  +  uf+1 )  -  (uf+1  +  uj* )  _ 


2At 


e  i  (^”+Y  -  -of+} ) + (1  -  9,  )(uf+1  -  of) 


- r  L 

'e,( 

— 

n+l  —  n+l') 

Hj+l  -Tlj  J 

Ax 

+  (l-0i) 

(<+!-<) 

Ax 

+2 


e«-dj"+1)+(i-e)(df+1-df)' 

Ax 


-gs  o 


y(l-S;)/c0AJCp(l-p)r| _ ^ 

B  8^ 


(u-u)2  =  0. 


(129) 
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Solution  of  the  system  of  equations 

A  Newton-Raphson  iteration  procedure  was  used  to  solve  the  system  of  equa¬ 
tions  given  above.  If  an  equation  is  designated  as 

F  =  F(dj ,  T|j ,  Uj ,  \>j ,  dj+i ,  r|j+i ,  Wj+| ,  ^j+i  j  (130) 

then  for  any  two  points  'j'  and  'j  +  Y  at  any  time,  it  can  be  expanded  as  a  Taylor 
Series 


m+lp_mp  +  dF  ■ 

dd ; 


9F 


j+i 


Aij+1  + 


9F 

dll: 


AU>  +  dll: 


dF 


1+1 


Au, 


1+1' 


dF 

90; 


dF 
dr\j 

a  dF 

Av>+*r 


ATlj 


dF 


dr\j+: 


Arjj+i  ■ 


j+i 


Au 


(131) 


l+i 


where  'm'  and  'm  +  V  indicate  the  iteration  level  and  'j'  and  'j  +  V  indicate  the  x- 
location.  The  partial  derivatives  are  evaluated  on  the  basis  of  the  values  of  the  vari¬ 
ables  after  the  'mth'  iteration.  Therefore,  each  equation  can  be  transformed  into  a 
linear  equation  in  terms  of  the  eight  unknowns:  Ady  Ary,  Aty,  Ary ,  Adj+1,  Ary+1,  A iy+1, 
and  Aa)j+1.  The  goal  of  the  method  is  to  solve  for  the  change  between  'm'  and  'm  +  V 
such  that m  +  2F  =>  0.  It  can  also  be  seen  that  mF  is  a  function  of  the  final  values  of  the 
variables  at  time  'ri  and  the  values  for  the  'mth/  iteration  at  time  'n  +  1/  which  are 
all  known  values.  The  transformed  equations  are  then  put  into  the  form 


=  flArfj  +  bAt|j  +  cAuj  +  dAUj  +  eAij+1  +/Ar)j+1  +  gAUj+1  +  FAr>j+1 


where 


(132) 


a  = 


e  = 
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dF 


dF 
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UI 

C=  9m”  / 

dOj  ' 

dF 

9F 

dF 

T  ^F 

ddi+ 1  ' 

9hj+i  ' 

g  =  a 

0Mj+1  ' 

h=  a 
9nj+i 

,m+ip_mp/  and  the  goal  is  to  achieve  w+l p  q 

K  =-  F 


(133) 


For  example,  the  conservation  of  water  mass  equation  becomes 


2A  t 


e{dft-d^)+(i-Q)(d»+1-d?) 


Ax 


e(wf+V  -  uf+1) + (i  -  e)(uf+1  -  ) 


Ax 


=  0 


(134) 


and  the  derivatives  are  evaluated  as 
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_0F__1 _ ud  0 

ddj  2At  Ax  2 


^-Mf+1)  +  (1-0)(»f+1-»f) 


3F  _0 

dlti  2 


i-f-o 

oDj 


3F  _  i  ,  ue  (  e 

dd]+i  2A  t  Ax  2 


/=^L-=° 

^ij+i 


3F  e  r nil] ,  de 

dwi4.i  2  Ax  Ax 


.  dF  n 

/Z  =  ^“  =  0-  (142) 

The  same  procedure  was  followed  for  the  other  three  equations,  using  the  nota¬ 
tions  F//F",F,"/a'/a'/,a'"/.../K,,K"/K'"  to  distinguish  among  equations.  The 
discretized  equations  and  their  derivatives  are  listed  in  Appendix  A. 

For  each  point,  there  are  now  four  unknowns  {Ad,  Arj,  A u,  and  An)  or  a  total  of 
4 N  unknowns,  where  N  is  the  number  of  cross  sections.  For  each  reach  there  are 
four  linear  equations  or  4 (N  - 1)  equations.  Four  boundary  condition  equations  are 
needed  to  close  the  system.  The  linear  algebraic  system  can  be  represented  as 


[A]{AZ}  =  {K]  or  {AZ}  =  [A]-1{K} 


where 


{AZ}  =  {Ad1 ,  Arij ,  Aux ,  Ao^ ,  Ad2 ,  Ar|2 , .. AdN ,  At|n  ,  Amn  ,  AoN 


{K}  =  {KlrK[,  K'{,  K’{’  K2,K2,  K2,  K2,  ...,Kn,KU,K'^,  K$}T .  (145) 

The  coefficient  matrix  [A]  is  depicted  in  Figure  31.  [A]  is  a  banded  coefficient 
matrix  of  11  diagonals  with  the  x's  signifying  boundary  conditions  and  the  empty 
spaces  filled  with  zeros.  The  system  is  solved  using  a  decomposition-back- 
substitution  scheme  to  obtain  the  A 2  terms.  Subsequently,  for  the  next  iteration 
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Figure  31.  Coefficient  matrix  for  the  Newton-Raphson  iteration  procedure. 


m+1d?+l=md?+l+Ady  (146) 

The  values  of  the  variables  are  updated  in  this  way  for  the 'm  +  Y  iteration  and 
solved  again  until  a  specified  tolerance  is  met.  The  tolerance  is  expressed  in  a  least 
squares  form 

N- 1/  -  \2  N  ~ 

Tol=  I  m+1F-mF  =  I(A's)2.  (147) 

When  the  tolerance  is  met,  the  final  'm  +  V  iteration  values  are  used  as  the  values  of 
the  dependent  variables  for  time  'n  +  Y.  The  time  is  then  incremented  and  the 
iteration  procedure  begins  again.  For  the  first  iteration  at  a  time  step  (' m  =  O'),  the 
initial  values  of  the  variables  are  taken  to  be  equal  to  the  final  value  of  the  variables 
at  the  previous  time  step. 

Ice  cover  stability,  solution  methods,  and  boundary  conditions 

The  conservation  of  ice  momentum,  as  expressed  in  eq  102  or  115,  is  actually  an 
equation  expressing  a  unique  balance  of  ice  momentum  and  external  forces  acting 
on  an  ice  jam  element.  Through  assumption  of  the  passive-pressure  failure  criteria 
and  adopting  the  Rankine  equation  for  (the  passive  pressure  coefficient),  the 
equation  is  valid  for  a  static  element  of  an  ice  jam  at  its  limit  of  stability  The  equa¬ 
tion  is  also  valid  for  an  ice  jam  element  failing  and  in  motion.  In  the  former  case,  if 
the  forces  on  a  static  element  exceed  the  passive  pressure  limit,  application  of  the 
equation  (in  concert  with  the  other  conservation  equations)  results  in  downstream 
ice  movement.  For  an  ice  element  in  motion,  application  of  the  equation  for  a  change 
in  the  forces  may  result  in  changes  in  ice  velocity  or  thickness.  For  many  cases,  the 
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net  force  applied  to  the  ice  cover  element  may  be  less  than  that  "allowed"  by  the 
passive  pressure  limit.  Application  of  the  equation  in  this  case  would  result  in 
either  ice  motion  upstream  (negative  ice  momentum)  or  thinning  of  the  jam  to 
achieve  a  balance.  In  natural  systems,  thinning  or  upstream  movement  of  ice  rarely 
occurs.  The  forces  exerted  on  a  thickened,  but  static,  section  of  jam  are  reduced. 
The  jam  simply  remains  in  place  with  no  change  in  thickness.  In  such  a  case  when 
the  forces  exerted  lead  to  a  stable  jam  (i.e.,  below  the  passive  pressure  limit),  the 
conservation  of  ice  momentum  equation  must  be  adjusted  so  that  negative  ice 
velocities  or  thinning  do  not  take  place. 

The  computational  procedure  commences  at  the  upstream  end  of  the  system, 
determining  the  stability  of  each  reach  at  the  beginning  of  each  time  step.  The  forces 
acting  on  the  jam  and  water  in  the  reach  are  assessed  and  compared  with  the  maxi¬ 
mum  resistive  force  at  the  downstream  end  of  the  reach.  The  assumption  of  passive 
pressure  limits  the  net  force  at  the  downstream  end  of  each  reach  to  that  given  by 
eq  89.  If  the  sum  of  the  forces  (normal  force  at  the  upstream  end  of  the  reach,  ice 
momentum,  hydrostatic  pressure,  gravity,  shear  on  the  underside  of  the  cover,  and 
friction  at  the  banks)  is  less  than  the  passive  pressure  limit,  the  reach  is  deemed 
stable.  The  stability  check  continues  in  the  downstream  direction  with  the  normal 
force  at  the  upstream  end  of  the  next  reach  set  equal  to  the  net  force  at  the  down¬ 
stream  end  of  the  current  reach.  The  stability  check  assesses  the  stability  of  each 
cross  section  by  comparing  the  net  force  and  passive  pressure  resistance. 

Modifications  to  the  conservation  equations  of  ice  mass  and  ice  momentum 
applied  to  each  reach  are  necessary,  depending  on  the  stability  conditions  at  the 
ends  of  each  reach.  One  of  the  upstream  boundary  conditions  of  the  system  pro¬ 
vides  a  relationship  for  the  upstream  ice  thickness  of  the  first  reach.  Thus,  the  con¬ 
servation  of  ice  mass  equation  provides  a  relationship  for  the  upstream  ice  velocity 
of  a  reach.  The  conservation  of  ice  momentum  provides  a  relationship  for  the  down¬ 
stream  ice  thickness  for  a  reach.  One  of  the  downstream  boundary  conditions  of 
the  system  describes  a  relationship  for  the  downstream  ice  velocity  at  the  last  reach. 
If  a  cross  section  is  deemed  to  be  stable,  ice  velocity  in  it  will  be  zero,  and  the 
conservation  of  ice  mass  equation  reverts  to 

uf+1=0.  (148) 

Several  combinations  of  upstream  and  downstream  stability  are  possible  for  a 
reach  (Fig.  32).  For  the  first  case,  both  upstream  and  downstream  ends  are  unstable 
and  the  full  equations  are  used.  In  the  next  case,  both  ends  are  deemed  stable  and, 
thus,  the  ice  velocity  equals  zero  at  both  ends,  and  eq  148  is  used.  Because  there  is 
no  ice  movement,  the  ice  thickness  must  remain  constant.  Consequently,  the  con¬ 
servation  of  ice  momentum  equation  reverts  to 

<+Y-Tif+1=0.  (149) 

In  the  third  case  of  Figure  32,  the  downstream  end  is  stable  and,  thus,  ice  velocity  is 
equal  to  zero.  But  the  upstream  end  is  unstable  and  the  full  ice-continuity  equation 
is  required.  Ice  moves  into  the  reach,  thickening  at  the  downstream  end,  and 
requiring  use  of  the  full  ice  momentum  equation.  The  last  case  shown  occurs  when 
the  ice  is  stable  at  the  upstream  end,  but  unstable  at  the  downstream  end.  Again,  eq 
148  is  used  for  the  ice  continuity  equation,  but  the  full  ice  momentum  equation  is 
required  for  defining  the  downstream  thickness. 
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Figure  32.  Conditions  of  ice  jam  stability  for  a  computational  reach. 


Three  modes  of  solution  of  the  system  of  equations  are  presented  below.  The 
fully  coupled  mode  comprises  the  simultaneous  solution  of  all  four  dependent 
variables  using  the  full  equations  and  Newton-Raphson  iteration  technique  pre¬ 
sented  above.  The  uncoupled  or  loosely  coupled  mode  solves  the  water  equations 
and  ice  equations  separately  and  sequentially,  as  is  done  in  the  moving  ice  model 
of  Tsai  et  al.  (1988).  The  static-unsteady  thickness  mode  first  solves  the  water  equa¬ 
tions,  then  a  modified  version  of  the  conservation  of  momentum  equation  not 
including  the  ice-momentum  component.  This  latter  procedure  is  comparable  to 
the  procedure  used  in  stationary  jam  formulations. 

Fully  coupled  solution 

The  fully  coupled  mode  of  solution  uses  the  full  equations  of  mass  and  momen¬ 
tum  conservation  for  the  water  and  ice  with  the  simplifications  mentioned  above 
for  sections  and  reaches  that  are  determined  to  be  stable.  The  four  dependent  vari¬ 
ables  are  solved  simultaneously  using  the  Newton-Raphson  iteration  scheme  men¬ 
tioned  above.  The  coefficient  matrix,  as  depicted  in  Figure  31,  is  solved  by  decom¬ 
position-back-substitution  using  banded-matrix  manipulation  techniques.  The 
coefficient  matrix  is  a  banded  symmetrical  matrix  with  five  bands  above  and  five 
below  the  main  diagonal.  Four  equations  for  each  reach  relate  the  4 N  unknowns, 
where  N  is  the  number  of  cross  sections.  Two  boundary  condition  equations  at  the 
upstream  end  and  two  at  the  downstream  end  are  required  to  close  the  system.  The 
upstream  boundary  conditions  are  typically  a  water  discharge  relation  (e.g.,  speci¬ 
fied  discharge)  and  an  ice  thickness  relation  (such  as  static  equilibrium  or  specified 
thickness).  The  downstream  boundary  conditions  are  typically  a  specified  depth 
relation  (under  ice)  and  an  ice  velocity  relation. 
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Figure  33.  Block  diagram  for  the  fully  coupled  solution  scheme. 


The  block  diagram  for  the  fully  coupled  solution  scheme  is  presented  in  Figure 
33.  The  initial  conditions  are  read  and  the  timer  incremented.  Boundary  conditions 
for  the  new  time  step  are  read  and  local  jam  stability  checked.  The  coefficient 
matrix  is  then  computed  and  solved,  tolerances  are  checked,  and  variables  are  reset 
for  the  next  time  step.  For  special  cases  when  the  entire  ice  cover  is  deemed  stable, 
the  solution  reverts  to  only  the  water  variables,  as  the  ice  variables  will  remain 
constant.  In  this  case,  an  abbreviated  coefficient  matrix  is  developed  and  solved 
that  contains  only  two  equations  per  reach  and  the  two  water  boundary  condition 
equations  (water  discharge  upstream  and  water  depth  downstream).  This  abbrevi¬ 
ated  matrix  is  a  banded  symmetrical  matrix  with  two  bands  above  and  two  below 
the  main  diagonal.  The  reversion  to  the  abbreviated  coefficient  matrix  reduces  com¬ 
putation  time. 

Loosely  coupled  solution 

The  loosely  coupled  mode  of  solution  also  uses  the  full  equations  of  mass  and 
momentum  conservation  for  the  water  and  ice,  including  the  simplifications  for 
sections  and  reaches  that  are  determined  to  be  stable.  However,  in  this  mode,  the 
water  variables  are  solved  separately  and  sequentially  from  the  ice  variables.  Both 
solutions  are  by  the  Newton-Raphson  iteration  scheme  presented  above.  The  coef¬ 
ficient  matrix  for  each  solution  is  a  banded  symmetrical  matrix  with  two  bands 
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Figure  34.  Block  diagram  for  the 
loosely  coupled  solution  scheme. 


above  and  two  below  the  main  diagonal.  Two  equations  for  each  reach  relate  the 
2 N  (water  or  ice)  unknowns.  A  single  boundary  condition  at  the  upstream  and  the 
downstream  ends  is  required  to  close  the  system. 

The  block  diagram  for  the  loosely  coupled  solution  scheme  is  presented  in  Fig¬ 
ure  34.  The  initial  conditions  are  read,  the  timer  incremented,  and  boundary  condi¬ 
tions  for  the  new  time  step  obtained.  The  ice  variable  solution  can  have  a  different 
time  step  than  the  water  variable  solution  and  is  specified  as  an  even  fraction  of  the 
water  time  step.  This  procedure  is  needed  because  the  ice  variables  can  change 
faster  than  the  water  variables,  especially  when  moving  ice  stops  and  thickens.  The 
ice  cover  stability  is  checked,  the  ice  variable  coefficient  matrix  is  computed  and 
solved,  tolerances  are  checked,  and  the  ice  variables  are  reset  for  the  next  ice  time 
step.  If  the  entire  ice  cover  is  found  stable,  the  ice  variables  are  simply  reset  and  the 
ice  time  step  incremented.  Following  the  correct  number  of  ice  time  step  calcula¬ 
tions,  the  water  coefficient  matrix  is  computed  and  solved,  tolerances  checked,  and 
the  water  variables  reset.  This  solution  mode  is  included  as  a  comparison  to  the 
fully  coupled  mode  because  there  are  currently  no  existing  formulations  that  are 
fully  coupled. 
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Static-unsteady  thickness  solution 

This  solution  mode  is  the  same  technique  used  by  several  models  where  the 
unsteady  water  flow  is  solved,  followed  by  a  calculation  of  the  ice  thickness.  In  this 
mode,  the  water  coefficient  matrix  is  solved  and  then  the  ice  thickness  calculated 
by  an  abbreviated  form  of  the  conservation  of  ice  momentum  equation,  in  which 
ice  velocity  is  always  assumed  to  be  zero.  The  resulting  calculation  of  ice  thickness 
may  violate  ice  mass  conservation,  since  the  ice  velocity  is  neglected  and  is  an 
instantaneous  (not  time-dependent)  reassessment  of  ice  thickness.  Equation  129  is 
modified  to 


(g(l-Si)Kp(l-p)  +  gSi) 


Hff-nr1) 

{  rr 

Ax 

+s 

Ax 

-£S0 


g(l-Si)Vt£p(l-p)Tl _ /i_^  =  0 

B  8Sj  fj  “ 


(150) 


where  the  average  (bar)  terms  are  given  as  simple  x-dependent  averages.  Previous 
formulations  for  jam  thickness  by  a  relation  similar  in  form  to  eq  150  indicate  that 
the  integration  can  move  either  in  an  upstream  or  downstream  direction  (Beltaos 
1993).  This  study  found,  however,  well  behaved  solutions  only  for  integration  in 
the  downstream  direction,  given  a  thickness  relation  for  the  upstream  boundary. 
An  equilibrium  thickness  condition  is  specified  at  the  upstream  boundary  by  set¬ 
ting  the  d/dx  terms  of  eq  150  to  zero  and  solving  the  resulting  relation  for  T|.  In 
recognizing  that  the  thickness  does  not  decrease  upon  a  reduction  in  the  forces 
acting  on  the  cover,  the  newly  calculated  thickness  is  compared  to  the  previous 
thickness  and  the  greater  value  adopted  as  the  new  thickness.  This  solution  mode 
is  the  least  computationally  intensive  but  neglects  the  effects  of  ice  momentum.  It 
is  included  in  this  study  only  for  comparison  to  the  loosely  coupled  and  fully  coupled 
modes. 


Boundary  conditions 

Boundary  condition  equations  are  necessary  for  the  closure  of  the  system  of  equa¬ 
tions  for  each  solution  technique  described  above.  In  general  terms,  upstream  and 
downstream  boundary  conditions  are  specified  in  each  case.  For  the  fully  coupled 
solution  mode,  ice  and  water  boundary  conditions  are  necessary  at  the  upstream 
and  downstream  boundaries.  Similar  to  open-water  modeling,  the  water  bound¬ 
ary  conditions  are  typically  a  specified  water  discharge  at  the  upstream  and  a  speci¬ 
fied  depth  relation  at  the  downstream  boundaries.  The  ice  boundary  conditions 
include  ice  thickness  at  the  upstream  boundary  and  ice  velocity  at  the  downstream 
boundary.  As  mentioned  in  the  foregoing  sections  describing  the  various  solu¬ 
tion  techniques,  various  relations  can  be  used  to  describe  these  specified  boundary 
conditions. 
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NUMERICAL  MODEL  DESCRIPTION 


Model  description  and  capabilities 

A  numerical  model,  called  THKFUL,  was  developed  on  the  basis  of  the  one¬ 
dimensional,  fully  coupled  equations  presented  in  the  Formulation  section.  The 
model,  written  in  the  FORTRAN  programming  language,  computes  the  unsteady 
water  depth,  velocity,  ice  thickness,  and  ice  velocity  for  a  river  channel.  The  current 
version  of  THKFUL  assumes  a  prismatic,  rectangular  channel  with  uniform  bed 
slope,  constant  values  of  bed  and  ice  roughness,  and  invariant  and  uniform  ice 
properties. 

The  model  THKFUL  calls  on  many  subprograms  to  do  various  computational 
tasks.  A  file  utility  subprogram  opens  all  required  input  and  output  files.  Another 
subprogram  reads  all  the  channel  geometry,  ice  property,  and  physical  constant 
data.  One  subprogram  reads  the  initial  conditions  file  of  depth,  velocity,  ice  thick¬ 
ness,  and  ice  velocity  at  each  cross  section.  Another  reads  the  new  boundary  condi¬ 
tion  data  at  the  beginning  of  each  time  step.  The  stability  of  the  ice  cover  at  each 
cross  section,  assessed  by  means  of  a  force  accumulation,  is  determined  using 
another  subprogram.  One  of  two  equation-system  solvers  is  called  to  compute  flow 
and  jam  variables.  FULSOL  fills  the  coefficient  matrix  using  the  technique  described 
in  the  Formulation  section  for  the  simultaneous  solution  of  the  four  dependent  vari¬ 
ables.  However,  if  the  entire  ice  cover  is  found  to  be  stable,  only  the  water  variables 
are  calculated  by  WATUNC,  which  fills  and  solves  a  smaller  coefficient  matrix. 
Another  subprogram  writes  data  to  three  output  files  with  user-definable  write 
intervals.  The  first  output  file  provides  a  listing  of  the  depth,  velocity,  ice  thickness, 
and  ice  velocity  at  each  cross  section  at  the  specified  print  interval.  The  second 
output  file  provides  profile  data  for  the  bed,  bottom,  and  top  of  the  ice  cover,  and 
water  level  at  the  specified  print  interval.  The  third  output  file  contains  data  neces¬ 
sary  for  compiling  animated  plots  of  the  profiles  with  time.  A  full  program  listing 
can  be  found  in  Zufelt  and  Ettema  (1996). 

Another  program,  called  UNCTHK,  is  very  similar  to  THKFUL,  but  operates  in 
a  loosely  coupled  mode.  The  ice  and  water  variables  are  calculated  separately  and 
sequentially  In  this  model,  the  ice  variable  solver,  ICEUNC,  is  called  one  or  more 
times  (as  specified  by  the  user)  for  each  time  that  the  water  solver,  WATUNC,  is 
called.  This  modification  was  made  because  the  ice  variables  can  change  rather 
abruptly,  e.g.,  when  a  moving  ice  cover  stops  and  thickens.  The  water  variables 
tend  to  respond  slower  and  smoother.  One  would  expect  that  the  results  obtained 
from  the  loosely  coupled  and  fully  coupled  models  would  approach  each  other  as 
the  time  step  is  reduced.  A  program  listing  for  UNCTHK  is  also  included  in  Zufelt 
and  Ettema  (1996). 

The  remainder  of  this  section  demonstrates  the  capabilities  of  the  model  THKFUL 
with  sample  output  plots  for  a  "baseline"  configuration  of  geometric,  hydraulic, 
and  ice  characteristics.  The  robustness  of  the  calculation  technique  is  demonstrated 
and  the  results  of  sensitivity  testing  with  relation  to  Courant  number  and  theta 
weighting  factors  for  the  water  and  ice  are  presented.  Alternate  boundary  condi¬ 
tions  are  described  and  the  effects  of  using  variable  size  length  steps  are  addressed. 

Baseline  runs 

A  baseline  configuration  was  developed  to  provide  a  standard  against  which 
future  runs  could  be  compared.  Care  was  taken  in  developing  this  5-km-long 
"hypothetical  river"  to  make  sure  that  the  geometric,  hydraulic,  and  ice  character- 
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istics  were  in  the  middle  of  the 
ranges  expected  for  ice-covered 
rivers  at  the  time  when  ice  jam¬ 
ming  would  be  expected.  Table  1 
provides  the  values  of  several  key 
parameters  chosen  for  the 
baseline  configuration.  The  val¬ 
ues  of  the  Darcy-Weisbach  fric¬ 
tion  factor  were  chosen  to  repre¬ 
sent  a  system  where  the  ice  jam 
roughness  was  significantly 
greater  than  the  bed  roughness. 
Values  of  Kp  and  \i  were  obtained 
as  explained  in  the  Laboratory  Ex¬ 
periments  section.  Using  these  val¬ 
ues  and  eq  27  leads  to  kfk  =  0.325.  By  substitution,  the  values  in  Table  1  for  these 
two  parameters  were  chosen. 

In  addition  to  the  parameters  in  Table  1,  the  tolerance  for  each  variable,  as 
expressed  in  eq  147,  was  set  to  lCL4  m.  The  upstream  boundary  conditions  were 
specified  as  water  discharge  and  equilibrium  jam  thickness  by  eq  147  with  the  non- 
uniform  terms  set  to  zero.  The  downstream  boundary  conditions  were  specified  as 
normal  depth  beneath  a  uniformly  thick  ice  cover  and  zero  ice  velocity.  The  up¬ 
stream  water  discharge  for  the  baseline  testing  begins  at  a  steady  value  of  100 
m3/s  for  20  minutes,  rises  to  200  m3/s  during  the  next  20  minutes,  and  remains  at 
that  level  for  a  total  test  time  of  600  minutes.  This  upstream  discharge  hydrograph 
is  shown  in  Figure  35.  The  maximum  number  of  iterations  of  the  Newton-Raphson 
solution  scheme  per  time  step  was  set  at  four.  However,  rarely  more  than  two  itera¬ 
tions  were  required  to  reach  the  specified  tolerance. 

The  initial  conditions  for  the  baseline  runs  were  determined  by  running  the  model 
with  a  steady  upstream  water  discharge  of  100  m3/s  and  a  prescribed  uniform 
thickness  jam  that  was  significantly  greater  than  the  equilibrium  thickness  pre¬ 
dicted  from  eq  25  to  ensure  that  there  would  be  no  further  thickening  as  the  water 
depths  and  velocities  adjusted  to  their  steady,  uniform  values.  Once  the  initial  con¬ 
ditions  for  water  depth  and  velocity  were  determined,  the  initial  uniform  ice  thick¬ 
ness  for  the  baseline  testing  was  chosen  to  be  only  slightly  greater  than  the  equilib¬ 
rium  value  calculated  using  eq  25.  For  the  baseline  rims,  this  procedure  resulted  in 
the  following  set  of  initial  uniform  values  of  water  depth  and  velocity,  ice  thick¬ 
ness,  and  ice  velocity:  1.729  m,  0.578 
m/s,  1.50  m  and  0  m/s,  respectively.  The 
equilibrium  jam  thickness  was  calcu¬ 
lated  to  be  1.47  m  from  eq  25. 

One  of  the  most  difficult  tasks  of  mod¬ 
eling  jam  shoving  and  thickening  is  the 
presentation  of  results.  The  process  is 
highly  unsteady,  with  large  variations  in 
the  values  of  the  dependent  variables 
over  fairly  short  periods.  Also,  consid¬ 
ering  the  interest  in  all  four  of  the  vari¬ 
ables  mentioned  above  means  that  a 
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Figure  35.  Upstream  water  discharge 
hydrograph  for  baseline  run. 


Table  1.  List  of  baseline  testing  parameters. 


Number  of  cross  sections 

101 

Length  step 

100  m 

Time  step 

30  s 

Channel  width 

100  m 

Bed  slope 

0.0005 

Solid  ice  specific  gravity 

0.92 

Ice  cover  porosity 

0.40 

Darcy-Weisbach  friction  factor — bed 

0.08 

Darcy-Weisbach  friction  factor — ice  cover 

0.12 

kQ — coefficient  of  lateral  pressure 

0.50 

X — coefficient  of  friction,  ice-on-ice 

0.65 

Kp — passive  pressure  coefficient 

3.85 

Theta  weighting  factor — water  (0) 

0.60 

Theta  weighting  factor — ice  (0j) 

0.60 

Maximum  number  of  iterations /time  step 

4 
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single  plot  of  them  becomes  overly  cumbersome  to  interpret.  The  interaction 
of  the  dependent  variables  with  and  on  each  other  necessitates  a  form  of  output 
illuminating  the  changes  in  all  the  variables  with  time.  One  possible  format 
is  to  present  the  output  in  the  form  of  profile  plots  at  specific  times.  While  these 
sequential  time  plots  can  be  animated  with  additional  graphic  enhancements,  they 
cannot  be  presented  easily  in  text  form.  Also,  as  the  length  of  the  modeled  system 
increases,  changes  of  values  in  the  vertical  become  less  obvious  because  of  loss  in 
resolution.  Therefore,  the  water  variables  (depth  and  velocity)  are  plotted  on  one 
set  of  axes  directly  above  the  ice  variables'  (thickness  and  velocity)  axes.  This  is 
done  for  several  times  at  points  on  the  inflow  hydrograph  where  there  are  signifi¬ 
cant  changes.  Figure  36  presents  the  results  of  the  baseline  model  for  several  time 
slices. 

Many  interesting  observations  can  be  made  from  Figure  36.  By  25  minutes  (5 
minutes  after  the  initial  increase  in  upstream  discharge),  the  water  velocity  has 
increased  near  the  upstream  end,  but  the  depth  has  not  changed  significantly.  The 
jam  is  thickening  somewhat  near  the  upstream  end,  but  the  ice  velocity  shows  the 
most  significant  change.  As  the  ice  moves,  the  shear  stress  between  the  ice  and  the 
water  decreases,  resulting  in  less  resistance  to  the  water  flow  and  lower  water  lev¬ 
els.  At  35  minutes,  ice  is  moving  throughout  the  system  and  is  by  no  means  uni¬ 
form.  The  ice  velocity  plot  shows  that  the  jam  failure  is  not  simply  a  progressive  or 
complete  failure  as  envisioned  in  the  laboratory  experiments,  but  is  rather  a  combi¬ 
nation  of  these  two  failure  modes,  with  several  areas  of  instability  arising  within 
the  jam.  Each  of  these  instabilities  moves  downstream  with  time,  resulting  in  local 
areas  of  thickening  where  ice  velocity  decreases.  The  stopping  and  starting  of  the 
ice  movement,  coupled  with  the  thickening,  results  in  changing  stress  levels  and 
local  states  of  stability  of  the  jam.  The  local  instabilities  move  downstream,  eventu¬ 
ally  dying  out.  At  45  minutes,  the  ice  at  the  upstream  end  of  the  system  has  ceased 
moving,  while  several  instabilities  are  continuing  to  move  downstream.  By  100 
minutes,  all  ice  motion  has  stopped  and  the  final  thickness  profile  prevails.  The 
water  depths  and  velocities  continue  to  adjust,  however,  until  they  reach  steady 
values. 

What  becomes  most  apparent  from  the  output  of  the  baseline  runs  is  the  high 
degree  of  interaction  between  the  dependent  variables.  There  is  no  simple  mecha¬ 
nism  by  which  an  ice  jam  fails,  moves,  and  thickens.  Shoving  and  thickening  are 
dominated  by  unsteadiness  and  nonuniformity.  Ice  velocity  depends  on  the  forces 
currently  exerted  against  the  jam  and,  thus,  is  changing  constantly  as  ice  comes 
into  motion  and  stops.  Local  changes  in  ice  velocity  or  thickness  cause  changes  in 
water  shear  on  the  underside  of  the  jam,  which  then  affect  the  forces  on  the  jam. 
The  final  thickness  profile  also  points  to  the  importance  of  ice  momentum.  The 
downstream  end  of  the  simulated  jam  had  a  boundary  condition  of  zero  ice  veloc¬ 
ity,  i.e.,  no  ice  momentum  at  that  location.  Inspection  of  the  ice  momentum  equa¬ 
tion  for  the  last  reach,  however,  shows  that  ice  thickness  at  the  downstream  end  is 
affected  by  the  change  in  momentum  between  the  last  two  cross  sections.  Ice  thick¬ 
ness  is  least  at  the  downstream  end,  where  the  effects  of  ice  momentum  are  least, 
and  is  1.73  m.  This  thickness  is  greater  than  the  equilibrium  thickness  (r\  =  1.70  m) 
calculated  using  eq  25,  which  does  not  include  the  effects  of  ice  momentum.  The 
upper  reaches  of  the  jam  have  greater  ice  velocities  (as  evident  in  the  plots  at  30  and 
35  minutes).  The  effects  of  arresting  ice  momentum  are  clearly  reflected  by  the  much 
higher  levels  of  jam  thickness  in  the  upstream  reaches. 
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Model  rigor 

An  unsteady  model  should  be  able  to  operate  in  steady  flow  mode  without  prob¬ 
lems.  This  requirement  is  demonstrated  in  the  previous  section,  when  steady  flow 
values  were  obtained  for  the  initial  conditions  of  the  variables.  A  model  should 
also  be  robust,  being  able  to  compensate  for  initial  conditions  far  from  steady  or 
uniform.  The  model  was  tested  under  a  variety  of  "bad  input"  conditions  to  see  if 
it  could  adjust  to  steady  flow  conditions. 

The  first  tests  involved  initial  conditions  that  included  either  a  significant  depth 
increase  or  decrease.  The  initial  cover  thickness  was  specified  such  that  no  thicken¬ 
ing  would  take  place.  The  objective  was  to  see  if  the  model  could  return  to  steady 


conditions  of  depth  and  velocity  (initial 
baseline  conditions)  given  these  nonuni¬ 
form  initial  conditions  of  depth.  The  up¬ 
stream  discharge  boundary  condition 
was  held  at  100  m3/s.  Figure  37  shows 
the  water  depth  and  velocity  at  0,  5,  10, 
20,  and  120  minutes  for  an  initial  local 
increase  in  depth.  By  5  minutes,  the  ini¬ 
tial  water  storage  "lump"  is  flattened  by 
gravitational  forces  acting  in  the  up¬ 
stream  and  downstream  direction.  It  then 
slowly  "drains"  through  the  downstream 
end  of  the  system  with  time,  reaching 
steady  uniform  flow  conditions  by  400 
minutes.  The  water  velocity  depicted  at 
time  0  is  that  associated  with  the  onset  of 
the  simulation  and  not  the  initial  condi¬ 
tions.  The  ice  thickness  remained  constant 
at  the  initial  value,  and  hence  ice  velocity 
remained  equal  to  zero  throughout  the 
run.  Figure  38  gives  the  initial  and  final 
water  surface  level  profiles  for  this  run, 
showing  the  impossible  nature  of  the  ini¬ 
tial  water  levels.  A  second,  similar  run 
looked  at  a  local  decrease  in  depth.  Fig¬ 
ure  39  shows  the  water  depth  and  veloc¬ 
ity  at  0, 5, 10, 20,  and  120  minutes.  In  this 
run,  the  initial  water  storage  deficit  sends 
negative  gravity  waves  both  in  the  up¬ 
stream  and  downstream  directions.  The 
deficit  slowly  fills  with  time,  reaching 
steady,  uniform  flow  conditions  by  400 
minutes.  Figure  40  shows  the  initial  and 
final  water-surface  level  profiles  for  this 
run. 

Another  run  was  made  using  initial  ice 
jam  thicknesses  that  would  result  in  im¬ 
possible  water  levels.  For  this  run,  the 
steady,  uniform  depth  of  the  baseline  con¬ 
dition  (1.73  m  at  a  steady  flow  rate  of  100 


Figure  37.  Water  depth  and  velocity  pro¬ 
files  at  various  times  for  the  initial  condi¬ 
tion  of  local  depth  increase. 


Figure  38.  Initial  and  final  water  surface 
level  profiles  for  depth  increase. 
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m3/ s)  was  overlaid  with  an  ice  jam  that 
was  generally  2.5  m  thick,  except  for  two 
areas  where  the  thickness  increased  (to 
3.0  and  2.8  m).  The  2.5-m  thickness  was 
chosen  so  that  the  jam  would  remain  sta¬ 
tionary.  The  water  depth  and  velocity 
beneath  the  jam  adjusted  to  its  final 
steady  values  by  150  minutes.  Figure  41 
shows  the  bottom  of  ice  and  water  sur¬ 
face  profiles  for  the  initial  and  final  con¬ 
ditions. 

The  final  check  on  model  rigor  used 
a  steady  water  flow  rate  of  100  m3/s,  but 
an  initial  ice  jam  thickness  that  was  less 
than  the  stationary  equilibrium  value  of 

l. 47  m  as  calculated  using  eq  25.  For 
these  runs,  the  ice  jam  was  initially  un¬ 
stable,  resulting  in  ice  movement  and 
subsequent  shoving  and  thickening  un¬ 
til  the  forces  on  the  ice  jam  attained  a 
balance.  Figure  42  shows  the  final  ice  jam 
thickness  profiles  for  runs  with  initial 
jam  thicknesses  of  1.3  and  1.45  m.  For 
the  initial  jam  thickness  of  1.45  m,  ice 
moved  only  in  the  upstream  reaches, 
with  the  ice  downstream  of  1500  m  re¬ 
maining  stationary  and  at  the  initial 
thickness.  The  shoving  and  resulting 
thickening  at  the  upstream  end,  where 
the  jam  grew  to  slightly  greater  than  1.47 

m,  caused  depth  to  increase  and,  thereby, 
shear  stress  on  the  ice  jam  underside  to 
decrease.  This  reduction  in  shear  stress 
ensured  the  stability  of  the  jam  in  the 
downstream  reaches. 

The  profile  for  the  initial  jam  thick¬ 
ness  of  1.3  m,  however,  shows  that  the 
jam  was  unstable  over  its  entire  length. 
Shoving  and  thickening  resulted  in  the 
final  jam  thickness  being  greater  than 


x  Location  (m) 

Figure  39.  Water  depth  and  velocity  pro¬ 
files  at  various  times  for  the  initial  condi¬ 
tion  of  local  depth  decrease. 


Figure  40.  Initial  and  final  water  surface 
level  profiles  for  depth  decrease. 


the  equilibrium  value  of  1.47  m  almost 
everywhere.  The  thinner  area  near  the  2000-m  location  is  caused  by  the  timing  of 
the  thickening  and  the  interaction  among  the  dependent  variables  in  this  fully 
coupled  solution.  The  stability  of  the  jam  at  any  section  depends  on  the  net  forces 
acting  on  the  jam,  which  in  turn  depend  on  the  levels  of  depth,  water  velocity,  jam 
thickness,  and  ice  velocity. 

This  check  confirms  the  potential  importance  of  ice  momentum  in  determining 
jam  thickness.  The  initial  jam  thickness  of  1.45  m  is  very  close  to  the  equilibrium 
value  but,  with  shoving  and  thickening,  the  eventual  thickness  is  slightly  greater. 
The  initial  jam  thickness  of  1.3  m,  while  fairly  close  to  the  equilibrium  jam  thick¬ 
ness  value,  results  in  a  rather  nonuniform  thickness  profile,  which  is  nearly  every- 
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where  greater  than  the  equilibrium  jam 
thickness  calculated  using  eq  25. 

Courant  number  sensitivity 

Many  calculation  schemes  are  sensi¬ 
tive  to  the  computational  time  and  space 
intervals  (A t,  Ax)  used.  The  sensitivities 
for  hyperbolic  problems  are  usually  dis¬ 
cussed  in  terms  of  the  Courant  number 


Figure  41.  Initial  and  final  water  surface 
level  and  bottom  of  ice  profiles  for  an  ice 
jam  having  nonuniform  thickness. 


Cr  = 


At\u  +  c\ 
Ax 


(151) 


where  c  is  the  gravity  wave  speed 
through  the  fluid.  Because  of  the  as¬ 
sumption  that  ice  floats  on  water  and  is 
free  to  move  up  and  down  in  the  verti¬ 
cal  direction  in  accordance  with  the  dic¬ 
tates  of  buoyancy  gravity  wave  speed 
beneath  a  jam  would  be  equivalent  to 
that  of  open  water,  i.e. 


■=4^- 


(152) 


Figure  42.  Final  jam  thickness  profiles  for 
two  initially  unstable  jam  thicknesses. 


The  Courant  number  expresses  the 
ratio  of  the  distance  traveled  by  a  dis¬ 
turbance  in  one  time  step  to  the  length 
of  a  computational  distance  step.  For  the 
simplest  Method  of  Characteristics,  the  Courant  number  must  be  less  than  or  equal 
to  unity  so  as  to  ensure  that  the  solution  remains  within  the  computational  do¬ 
main.  Implicit  finite  difference  methods,  as  well  as  interpolated-grid  Methods  of 
Characteristics,  relax  this  requirement  somewhat.  Care  must  be  taken,  however,  to 
make  sure  that  disturbances  do  not  travel  too  far  over  a  time  step,  thereby  resulting 
in  loss  of  resolution  or  too  much  smoothing.  For  a  given  computational  length  step, 
larger  Courant  numbers  imply  larger  time  steps  and  less  total  computation  time, 
with  concomitant  loss  of  temporal  resolution. 

Equations  151  and  152,  applied  to  the  baseline  parameters  in  Table  1  and  the 
initial  conditions  for  depth  and  water  velocity,  yield  Cr  =  2.82.  By  calling  this  value 
nominally  Cr  =  3,  runs  were  then  made  with  Cr  =  1,  5, 10, 12,  and  24.  The  Courant 
number  was  varied  in  these  tests  by  adjusting  the  time  step  and  leaving  the  length 
step  at  50  m.  The  boundary  conditions  file  was  modified,  however,  to  retain  the 
same  shape  and  timing  of  the  inflow  hydrograph.  Figure  43  shows  the  final  jam 
thickness  profile  for  the  runs  with  Cr  =  1,  3  (standard),  5, 10, 12,  and  24.  The  plot 
shows  that  there  is  little  difference  in  results  with  Cr  values  up  to  12,  and  that  some 
diffusion  or  smoothing  is  evident  for  Cr  =  24. 


Theta-weighting  factor  analysis 

The  theta-weighting  factors  for  the  water  0  and  for  the  ice  0t  are  presented  in  the 
Discretization  of  the  System  of  Equations  section  and  are  used  to  describe  the  time 
averaging  between  the  known  conditions  at  the  current  time  step  and  the  future. 
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Figure  43.  Final  jam  thickness  profiles  for 
Courant  numbers  of  1,  3,  5,  10,  12,  and 
24. 
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rather  suddenly.  The  effects  of  using  the  figure  44.  Final  jam  thickness  profilesfor  02 
current  values  of  ice  velocity  to  calculate  =  ^  an ^  ®  ~  ^5,  0.6,  0.66 ,  0.8,  and  1.0. 
the  new  values  were  unknown.  For  ex¬ 
ample,  if  the  ice  cover  was  in  motion  and  the  new  stability  check  showed  that  it 
should  be  stable  (the  new  ice  velocity,  X)  =  0),  the  previous  value  of  ice  velocity  may 
have  an  undesirable  effect  on  the  newly  calculated  ice  thickness  if  0*  <  1.0.  Detailed 
inspection  of  the  calculated  ice  velocities  (as  ice  begins  moving  and  as  it  slows  to  a 
stop),  however,  showed  that  absolute  changes  in  ice  velocity  over  a  time  step  were 
not  large  and,  thus,  the  0rweighted  average  of  ice  velocity  would  provide  a  proper 
solution. 

The  effects  of  the  water-weighting  factor  on  the  solution  were  investigated  by 
running  the  baseline  conditions  at  a  variety  of  0  values  for  0t  =  1.0.  Values  of  0 
tested  included  0.5, 0.55, 0.6,  0.66, 0.8,  and  1.0.  The  simulation  with  0  =  0.5  resulted 
in  unacceptable  oscillations  in  the  water  depth  and  velocity,  which,  because  the 
solution  is  fully  coupled,  resulted  in  large  oscillations  in  ice  velocity  and  thickness. 
The  oscillations  in  the  water  variables  aggravated  the  oscillations  in  the  ice  vari¬ 
ables  and  the  solution  became  highly  unstable.  For  instance,  the  depth  would 
decrease  with  an  accompanying  water  velocity  increase  that  increased  ice  velocity 
and  thickening,  thereby  further  reducing  flow  depth.  Figure  44  shows  the  final  jam 
thickness  profile  for  0  =  0.55,  0.6,  0.66,  0.8,  and  1.0.  While  0  =  0.55  shows  a  smooth 
final  profile,  it  resulted  in  highly  variable  water  depths,  especially  at  the  initial 
increase  in  upstream  water  discharge.  A  value  of  0  =  0.6  provided  the  most  accept¬ 
able  solution  in  terms  of  both  water  depth  and  ice  thickness  throughout  the  simu¬ 
lation  period. 
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Figure  44.  Final  jam  thickness  profilesfor  6j 
=  1.0  and  6  =  0.55,  0.6,  0.66,  0.8,  and  1.0. 
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The  effects  of  the  ice-weighting  fac¬ 
tor  were  then  investigated  by  running 
the  baseline  conditions  with  0  =  0.6  and 
0i  =  0.5,  0.55,  0.6,  0.66,  and  0.80.  Again, 
with  a  value  of  0*  =  0.5,  the  solution  be¬ 
came  highly  unstable  and  gave  unrea¬ 
sonable  results  in  terms  of  ice  thickness. 
Figure  45  shows  the  final  jam  thickness 
profile  for  0^  =  0.55, 0.6, 0.66, 0.8,  and  1.0. 
The  profile  shows  little  variation  with 
©i  =  0.6,  providing  physically  reasonable 
results.  Inspection  of  the  values  of  the 
other  dependent  variables  throughout 
the  simulations  showed  that  either  0j  = 
0.6  or  0.66  would  be  acceptable.  For  the 
remainder  of  the  model  testing,  the  val¬ 
ues  of  0  =  0.6  and  0j  =  0.6  were  used. 


Alternate  boundary  conditions 

The  boundary  conditions  for  the 
baseline  runs  included  specified  water 
discharge  and  equilibrium  thickness  at 
the  upstream  end,  and  zero  ice  velocity 
and  a  condition  of  normal  flow  depth 
beneath  the  jam  at  the  downstream  end. 
While  these  boundary  conditions  repre¬ 
sent  plausible  natural  conditions,  they 
are  by  no  means  all  inclusive.  A  dam,  for 
instance,  typically  has  a  water  surface 
slope  that  decreases  in  the  downstream  direction  and  results  in  downstream  water 
levels  that  are  significantly  above  the  normal  depth.  Ice  jams  in  dam  pools  often  are 
resistant  to  shoving  and  thickening  because  of  the  reduced  shear  stress  and  gravity 
forces  exerted  on  them.  The  equilibrium  jam  thickness  at  the  upstream  end  also 
implies  that  there  is  an  unlimited  supply  of  ice  upstream  of  the  jam  that  would 
continually  move  into  the  modeled  reaches  at  the  equilibrium  thickness. 

Two  alternative  boundary  condition  types  were  developed.  A  condition  of  speci¬ 
fied  ice  thickness  at  the  upstream  end  of  the  jam  facilitates  simulation  of  the  upper 
transition  zone  where  jam  thickness  is  reduced.  The  effects  of  lower  upstream  jam 
thickness  on  the  shape  of  the  final  jam  thickness  profile  can  also  be  investigated.  A 
condition  of  specified  water  depth  at  the  downstream  end  of  the  system  facilitates 
simulation  of  ice  jams  in  reservoirs  or  at  other  water-elevation  control  structures. 

A  run  was  carried  out  to  simulate  jam  shoving  and  thickening  in  a  reservoir 
where  the  downstream  depth  was  held  at  3.0  m.  The  initial  water  depths  and 
velocities  for  this  run  were  determined  by  running  a  steady  water  discharge  of  100 
m3/s  with  a  uniform  depth  throughout  the  3.0  m  and  letting  the  system  attain 
steady  flow  conditions.  The  baseline  inflow  hydrograph  was  then  run  with  all  other 
parameters  the  same  as  for  the  baseline  test.  Figure  46  shows  the  initial  and  final 
water  surface,  and  the  bottom  profiles  of  the  jams.  The  downstream  end  experi¬ 
enced  no  shoving  and  thickening,  remaining  at  the  initial  thickness  of  1.45  m,  while 


Figure  45.  Final  jam  thickness  profiles  for 
0=0.6  and  =  0.55, 0.6, 0.66, 0.8,  and  1.0. 


x  Location  (m) 


Figure  46.  Final  bed,  bottom  of  jam,  and  wa¬ 
ter  surface  level  profiles  for  the  condition  of 
downstream  depth  being  held  at  3.0  m. 
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the  upstream  reaches  thickened  to  over 
2.0  m.  The  initial  jam  thickness  profile  for 
this  test  was  specified  as  a  uniform  1.5  m. 

Given  the  water  surface  profile  that 
was  specified  as  the  initial  conditions,  the 
initial  stable  jam  thickness  would  have 
been  less  in  the  downstream  reaches.  This 
may  have  resulted  in  significant  shoving 
and  thickening  in  those  areas,  leading  to 
a  totally  different  final  jam  thickness  pro¬ 
file. 

For  another  run,  the  downstream 
depth  was  held  at  the  initial  uniform 
depth  of  1.729  m.  All  the  other  conditions 
were  the  same  as  the  baseline  run.  Fig¬ 
ure  47  shows  profiles  of  the  bed,  bottom 
of  jam,  and  water-surface  level  for  this 
run.  Confining  the  downstream  depth 
required  that  the  water  velocity  increase, 
which  resulted  in  significant  thickening 
and  water-surface  level  rises.  This  figure 
shows  that  the  profiles  of  the  water  sur¬ 
face  and  jam  bottom  at  the  downstream 
end  are  very  similar  to  the  classical 
downstream  transition  configuration  of 
an  idealized  ice  jam  (e.g.,  Ashton  1986). 

Several  other  runs  were  made  using  a 
specified  jam  thickness  as  the  upstream 
boundary  condition.  Three  runs  with  jam 
thicknesses  held  at  1.47, 1.4,  and  1.3  m  at 
the  upstream  end  were  compared  to  the 
baseline  run  (equilibrium  jam  thickness  at  the  upstream  end).  Figure  48  shows  the 
final  jam  thickness  profiles  for  these  runs.  The  profiles  are  very  similar,  with  differ¬ 
ences  primarily  occurring  in  the  most  upstream  1000  m.  The  implicit  nature  of  the 
numerical  solution  procedure  extends  the  effect  of  the  specified  (artificially  low) 
jam  thickness  for  several  computational  steps  downstream.  Also  evident  from  the 
figure  is  that  the  progression  of  the  shoving  and  thickening  differs  in  each  run, 
resulting  in  slightly  different  jam  thickness  profiles.  It  would  be  expected  that  a 
locally  weak  reach  would  shove  and  thicken  first,  setting  up  changes  in  the  water 
depth  and  velocity  (since  the  solution  is  fully  coupled).  Minor  changes  in  the  water 
variables  translate  into  minor  changes  in  the  ice  variables  and  in  the  final,  stable 
jam  thickness  profiles. 

Effects  of  variable  length  steps 

Final  testing  of  the  model  involved  examining  the  effects  on  the  predicted 
results  of  variable  computational  length  steps.  The  finite-difference  formulation 
allows  variable  length  steps,  whereas  other  computation  methods  do  not  without 
considerable  difficulties  (i.e.,  the  Method  of  Characteristics).  This  facility  is  needed 
because  field  data  are  often  variable  in  the  distance  between  cross  sections.  While 
the  computational  effects  of  variable  length  steps  should  be  minimal,  care  must  be 


Figure  47.  Final  bed ,  bottom  of  jam,  and 
water-surface  level  profiles  for  the  condi¬ 
tion  of  downstream  depth  being  held  at  an 
initial  depth  of  1.729  m. 


x  Location  (m) 

Figure  48.  Final  jam  thickness  profiles  for 
upstream  jam  thickness  being  held  at  1.47 , 
2.4,  and  1.3  m. 
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a.  In  the  first  1000  m  of  the  system . 


b .  In  the  last  1000  m  of  the  system. 


Figure  49.  Final  jam  thickness  profile  for  length  step  reduced  to  25  m. 


a.  Length  step  increased  to  100  m.  b.  Length  step  increased  to  200  m. 

Figure  50.  Final  jam  thickness  profile  for  length  step  changes  in  the  first  2000  m  of  the 
system. 


taken  not  to  violate  any  Courant  number  limitations  identified  earlier. 

The  baseline  run  had  uniform  length  steps  of  50  m,  which  give  Cr  =  3.  The 
Courant-number  sensitivity  findings  described  earlier  indicate  that  Cr  values  in 
the  range  of  1  to  12  should  provide  acceptable  results.  From  eq  151,  this  finding 
translates  to  acceptable  length  steps  ranging  from  16.7  to  200  m.  Several  runs  were 
made  with  a  variety  of  length  steps. 

The  first  run  involved  changing  the  length  step  to  25  m  for  the  first  1000  m  of  the 
system.  This  step  length  decreased  the  Cr  to  1.5  in  the  first  1000  m  of  the  system 
and  increased  the  total  number  of  cross  sections  to  121.  All  the  other  conditions 
were  the  same  as  for  the  baseline  run.  A  similar  run  decreased  the  length  step  to  25 
m  for  the  most  downstream  1000  m  of  the  system.  These  two  reaches  are  where  the 
greatest  variations  in  the  final  jam  thickness  profile  exist.  Figure  49  presents  the 
final  jam  thickness  profiles  for  these  two  runs,  compared  to  that  of  the  baseline  run. 
The  25-m  upstream  spacing  resulted  in  a  more  gradual  rise  to  the  maximum  thick¬ 
ness,  which  occurred  around  1000  m.  The  final  jam  thickness  was  up  to  0.05  m  less 
than  the  baseline  run,  but  only  for  locations  upstream  of  2000  m.  Almost  no  differ¬ 
ence  existed  between  the  baseline  and  the  25-m  downstream  spacing  runs,  espe¬ 
cially  in  the  most  downstream  1000  m  where  the  spacing  was  changed. 

The  next  two  runs  increased  the  length  step  in  the  most  upstream  2000  m  of  the 
system  to  100  and  200  m  respectively.  The  values  of  Cr  for  the  upstream  2000  m 
increased  to  6  and  12,  with  the  total  number  of  cross  sections  reduced  to  81  and  71 
respectively.  Figure  50  presents  the  final  jam  thickness  profiles  for  these  two  tests. 
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plotted  with  the  jam  thickness  profile  for 
the  baseline  test.  The  100-m  upstream 
spacing  resulted  in  a  final  thickness  pro¬ 
file  that  was  almost  identical  to  the 
baseline  run.  The  profile  shape  and 
maximum  thickness  were  very  similar. 
The  200-m  upstream  spacing  showed 
slightly  more  variation  than  the  100-m 
upstream  spacing  but  was  also  very 
similar  in  shape  and  maximum  thickness 
to  the  baseline  run. 

The  final  run  simulated  random  spac¬ 
ing  of  cross  sections  (typical  of  field  data) 
by  using  randomly  assigned  length  steps 
between  25  and  100  m  (intervals  of  5  m)  for  a  total  flow  length  of  5000  m  with  81 
cross  sections.  Now,  Cr  ranged  from  1.5  to  6.  Figure  51  shows  the  final  jam  thick¬ 
ness  profile  for  this  run  along  with  the  baseline  results.  Minor  variations  existed  in 
the  final  thickness  profile  (up  to  0.03  m),  primarily  in  the  most  upstream  2000  m  of 
the  system.  Overall  profile  shape  and  maximum  thickness,  however,  were  very 
similar  to  the  baseline  run.  This  check  on  model  robustness  shows  that  variable 
length  steps  have  minor  negligible  effects  on  predicted  results,  as  long  as  the  Cou- 
rant  number  limitations  identified  earlier  are  followed. 

Summary 

The  main  points  emerging  here  are  that  the  fully  coupled  model  of  ice  jam 
dynamics  is  robust  and  versatile.  Alternate  boundary  conditions  and  the  use  of 
variable  length  steps  make  the  model  adaptable  to  a  variety  of  physical  situations 
and  computational  capabilities.  The  sensitivities  of  results  to  computational  pa¬ 
rameters  are  found  to  be  minimal  as  long  as  the  parameters  are  held  within  reason¬ 
able  (practical,  physically  based)  ranges. 

Even  at  this  stage  of  use,  the  model  shows  that  the  effects  of  ice  momentum  on 
the  solution  are  quite  significant,  as  evidenced  in  the  profiles  of  jam  thickness.  For 
example,  the  baseline  jam  thickness  profile  shows  that  almost  the  entire  system 
ends  up  with  jam  thicknesses  greater  than  those  calculated  using  the  stationary 
jam  theory  embodied  in  eq  25.  Furthermore,  the  model  shows  that  shoving  and 
thickening  make  up  a  highly  unsteady  process,  marked  by  variations  and  interac¬ 
tions  between  the  dependent  variables  not  represented  in  past  modeling  efforts. 
Small  changes  in  any  of  the  dependent  variables  can  result  in  variations  in  the 
timing  of  the  shoving  and  thickening  event,  with  subsequent  variations  in  the  final 
jam  thickness  profile. 


Figure  51.  Final  jam  thickness  profile  for 
random  length  steps  throughout  the  system. 


UNSTEADY  JAM  DYNAMICS 
Effects  of  ice  momentum 

A  primary  objective  of  this  study  is  to  ascertain  the  effects  of  ice  momentum  on 
jam  formation  and  thickness.  The  importance  of  ice  momentum  already  became 
evident  when  the  stability  and  overall  robustness  of  the  numerical  model  devel¬ 
oped  for  the  study  were  evaluated.  The  fully  coupled,  numerical  model  includes 
ice  momentum  in  the  conservation  of  momentum  expression  formulated  for  ice  in 
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Figure  52.  Comparison  of  thickness  pro¬ 
file  predicted  by  static-unsteady  thickness 
model  to  fully  coupled  model  and  equilib¬ 
rium  thickness. 


a  jam.  The  static -unsteady  thickness 
solution  presented  earlier  allows  jam 
thickness  to  fluctuate  in  time  and  space, 
but  the  streamwise  movement  of  the  ice, 
and  therefore  ice  momentum,  are  not 
included. 

Early  formulations  of  equilibrium  jam 
thickness  are  based  on  uniform  flow 
depths  and  jam  thickness  in  a  reach  of 
constant  water  surface  slope.  As  a  result, 
the  force  balance  only  contains  the  fol¬ 
lowing  terms:  shear  stress  on  the  under¬ 
side  of  the  jam  from  the  water  flow,  the 
downstream  component  of  the  gravity 
force  on  the  jam  (based  on  the  bed  slope 
because  of  the  uniform  flow  assump¬ 


tion),  and  the  shear  resistance  at  the  banks.  Because  the  thickness  is  considered 
uniform  within  the  equilibrium  reach,  those  formulations  contain  no  terms  reflect¬ 
ing  changes  in  hydrostatic  pressure,  jam  strength,  water  depth,  or  jam  thickness 
along  the  reach  considered.  The  formulation  leading  to  eq  150,  however,  includes 
these  terms,  but  with  the  ice  velocity  terms  reduced  to  zero.  The  resulting  force 
balance  is  affected  by  the  nonuniform  thickness  of  the  jam  as  well  as  the  nonuni¬ 
form  depths  within  a  reach.  The  more  complete  balance  of  forces  expressed  in  eq 
150  enables  more  accurate  modeling  of  jam  profile,  whether  the  profile  includes  an 
equilibrium-thickness  section  or  not. 

To  show  the  effects  of  ice  momentum  on  jam  thickness  profile,  the  results  of  the 
static-unsteady  thickness  solution  leading  to  eq  150  were  compared  to  those 
obtained  by  the  fully  coupled  model.  The  loosely  coupled  moving-ice  model  was 
also  run  for  the  baseline  conditions  and  compared  to  the  fully  coupled  results  to 
demonstrate  the  effects  of  the  coupling  of  the  dependent  variables  on  the  jam  thick¬ 
ness  profile. 

The  static-unsteady  thickness  model  was  run  with  the  baseline  inflow  hydrograph 
and  all  other  parameters  the  same  as  for  the  baseline  condition.  Figure  52  shows 
the  predicted  jam  thickness  profile  obtained  by  the  static-unsteady  thickness  model. 
Also  plotted  are  the  thickness  profile  predicted  using  the  fully  coupled  model  and 
the  equilibrium  thickness  for  the  final  flow  level  of  200  m3/s.  The  static-unsteady 
force  balance  represented  by  eq  150  includes  the  partial  derivative  terms  represent¬ 
ing  changes  in  jam  strength  and  water  surface  slope  along  a  reach.  The  changing 
water  surface  slope  becomes  particularly  important  as  discharge  increases,  increas¬ 
ing  the  downstream-acting  forces  with  a  subsequent  increase  in  jam  thickness.  With 
a  milder  rate  of  increase  in  discharge,  the  water  surface  slope  would  remain  closer 
to  the  bed  slope,  resulting  in  a  static-unsteady  thickness  profile  nearer  to  the  equi¬ 
librium  thickness.  The  importance  of  ice  momentum  is  clearly  evident  from  the 
fully  coupled  thickness  profile  and  results  in  much  greater  ice  thickness.  Ice  veloc¬ 
ity,  and,  hence,  ice  momentum,  are  greatest  in  the  upstream  reaches  where  the 
increasing  discharge  wave  is  steepest.  The  ultimate  arrest  of  this  ice  momentum  is 
responsible  for  the  significantly  greater  thickness. 

The  effects  of  fully  coupling  the  solution  are  determined  by  comparison  with 
the  results  from  the  loosely  coupled  model.  The  loosely  coupled  moving-ice  model 
was  run  with  the  inflow  hydrograph  and  all  other  parameters  the  same  as  for  the 


baseline  condition.  The  number  of  ice 
calculation  cycles  for  each  water  calcu¬ 
lation  cycle  was  set  to  1.  Another  run  was 
made  with  the  number  of  ice  calculation 
cycles  set  to  2.  Figure  53  shows  the  final 
jam  thickness  profiles  for  the  fully 
coupled  and  loosely  coupled  model 
runs.  The  fully  coupled  model  results 
show  a  more  pronounced  effect  of  ice 
momentum.  As  the  number  of  ice  calcu- 


Figure  53.  Comparison  of  thickness  profile 
predicted  by  loosely  coupled  model  to  fully 
coupled  model 


lation  cycles  increases,  the  resulting 
thickness  profile  attains  a  shape  closer 
to  the  fully  coupled  thickness  profile. 

The  minimal  computational  time  saved 
by  using  the  loosely  coupled  model  with 

two  or  more  ice  calculation  cycles  is  outweighed  by  the  benefits  of  using  a  truly 
fully  coupled  model. 


Comparison  with  steady-state  models 

Steady-state  models  of  jam  thickness  existed  before  high-performance  mainframe 
and  personal  computers  became  available.  The  computational  power  that  exists 
today  allows  very  large  coefficient  matrices  to  be  solved  with  little  effort.  The 
execution  time  of  the  fully  coupled  model  for  tests  similar  to  the  baseline  configu¬ 
ration  described  in  the  Baseline  Runs  section  is  much  less  than  the  time  required  for 
compilation  of  input  files  and  analysis  and  plotting  of  output  files.  Even  though 
computational  speed  and  capability  have  greatly  increased,  there  still  are  reasons 
(such  as  unsteady  boundary  condition  information  requirements)  that  a  complex 
model  might  not  be  used.  For  many  situations,  a  steady-state  model  may  be  appro¬ 
priate  and  provide  results  within  a  specified  tolerance  or  accuracy  desired.  For 
instance,  Uzuner  and  Kennedy  (1976)  could  formulate  the  upstream  transition 
region  of  jam  thickness  using  a  frame  of  reference  that  traveled  upstream  at  the 
speed  of  the  progression  of  the  jam  front.  The  moving  reference  frame  made  the 
flow  quasi-steady  and  the  computation  much  easier. 

A  widely  used  steady-state  model  for  calculating  water  levels  in  ice-covered 
channels  is  HEC-2  modified  with  the  ice  cover  option  or  with  the  utility  ICETHK, 
or  both.  This  model  has  been  shown  (Zufelt  and  Doe  1986)  to  yield  accurate  results, 
provided  the  appropriate  values  of  the  jam-related  variables  are  chosen.  An  impor¬ 
tant  input  variable  for  a  steady-state  model  is  water  discharge.  However,  ice  jams 
are  markedly  unsteady  events  characterized  by  widely  fluctuating  discharges:  as 
level  ice  covers  (whence  the  jam  ice  originates)  and  ice  jams  fail,  releasing  stored 
water,  or  as  jams  form,  causing  local  reductions  in  discharge.  A  steady-state  model 
cannot  account  for  local  variations  in  discharge  caused  by  jam  formation  or  failure. 

The  discharge  record  shown  in  Figure  54  illustrates  how  unsteady  flows  can  be 
during  jam  breakup.  Fortunately,  the  stage  gauge  used  to  estimate  the  discharges 
in  Figure  54  was  located  approximately  1  km  downstream  of  the  toe  of  the  jam  and 
was  ice-free  during  the  event  (except  when  the  jam  failed  and  ice  passed  down¬ 
stream).  The  figure  clearly  shows  when  the  jam  initially  formed,  evidenced  by  a 
drop  in  the  discharge  downstream  at  the  gauge.  It  shows,  too,  the  rise  in  discharge 
as  the  water  level  in  the  jammed  reach  rose,  and  indicates  the  ultimate  failure  of  the 
jam,  evidenced  by  the  rapid  peak  and  passage  of  the  discharge  wave.  The  local 
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Figure  54.  Discharge  record  during  breakup  jam  initiation  and  failure. 


discharge  at  upstream  locations,  however,  could  be  markedly  different  from  that 
recorded  at  the  gauge.  In  many  cases,  gauges  are  not  located  near  the  jammed  reach 
and  estimates  of  discharge  become  difficult.  In  some  situations,  a  jam  may  create  a 
backwater  effect  on  the  measuring  gauge,  and  interpretation  of  the  stage  recording 
(in  developing  the  discharge  record)  also  becomes  questionable.  Winter  discharge 
estimates  can  be  off  by  as  much  as  25%  because  of  backwater  from  ice  jams.  As  a 
result,  the  choice  of  discharge  used  in  a  steady-state  model  can  greatly  affect  the 
results  of  the  modeling. 

The  fully  coupled  model  was  used,  in  conjunction  with  the  static-unsteady  thick¬ 
ness  model  with  a  constant  discharge,  to  determine  how  the  choice  of  water  dis¬ 
charge  affects  jam  thickness  profile.  Figure  55  shows  a  hypothetical  breakup  jam 
discharge  hydrograph.  The  hydrograph  is  characterized  by  a  fast  rise  and  very 
short  duration  peak,  similar  to  the  hydrograph  depicted  in  Figure  54.  While  the 
fully  coupled  model  uses  this  hydrograph  as  an  upstream  boundary  condition,  a 
steady-state  model  requires  that  a  single  discharge  be  chosen  as  an  input  variable. 
One  choice  might  be  the  peak  value  of  discharge.  If  the  modeled  system  were  of 
significant  length,  however,  the  instantaneous  peak  would  be  attenuated  as  it  trav¬ 
eled  downstream.  Thus,  the  discharge  felt  at  downstream  locations  (and  conse¬ 
quently  shear  stress  on  the  underside  of  the  jam)  would  be  less. 

Another  choice  for  the  steady-state  discharge  might  be  a  time-averaged  discharge 
based  on  the  time  required  for  a  disturbance  to  pass  through  the  system.  Two  tests 
were  run  with  the  static-unsteady  thickness  model  at  a  constant  discharge  level 
and  corresponding  uniform  values  of  depth  and  water  velocity.  This  simulated  the 
calculations  of  a  steady-state  model:  one  with  a  constant  discharge  of  the  peak 
value  of  200  m3/ s  from  Figure  55  and  the  other  with  a  time-averaged  constant 
discharge  value  of  167  m3/s. 

Figure  56  shows  the  final  jam  thickness  profiles  for  these  two  runs,  along  with 
the  final  jam  thickness  profile  predicted  using  the  fully  coupled  model  with  the 
unsteady  discharge  hydrograph.  The  steady  discharge  levels  and  uniform  initial 
conditions  result  in  uniform  jam  thickness  equal  to  that  given  by  the  equilibrium 
jam  formulation  represented  by  eq  25.  The  fully  coupled  thickness  profile  is  greater 
than  the  equilibrium  thickness  for  both  of  the  steady  discharge  choices  everywhere. 


71 


except  the  most  downstream  reaches. 
The  effects  of  ice  momentum  and 
nonuniformities  in  depth  and  thickness 
on  the  fully  coupled  profile  are  clearly 
evident,  even  for  this  discharge 
hydrograph  with  a  very  short  duration 
peak.  The  smaller  thicknesses  at  the 
downstream  end  are  ascribable  to  at¬ 
tenuation  of  the  peak  flow  as  it  travels 
downstream.  Figure  57  shows  the  dis¬ 
charge  hydrograph  at  the  upstream 
boundary,  the  midpoint  of  the  modeled 
system,  and  the  downstream  boundary 
The  peak  flow  is  greatly  reduced  as  it 
travels  downstream,  with  a  maximum 
value  less  than  150  m3/s  by  the  time  it 
reaches  the  midpoint  of  the  modeled 
system.  The  wave  is  also  significantly 
flattened.  The  high-frequency  fluctua¬ 
tions  in  discharge  at  the  mipoint  come 
from  ice  movement.  As  the  ice  moves 
and  stops,  resistance  to  water  flow  var¬ 
ies,  resulting  in  local  fluctuations  in  wa¬ 
ter  discharge. 
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Figure  55.  Hypothetical  breakup  jam  dis¬ 
charge  hydrograph . 
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Figure  56.  Final  jam  thickness  profiles  for 
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thickness  that  are  slightly  less  than  equi-  so  Ii'nInhIiiiiImiiIiihI  la-l±. 
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pected  flow  levels  (e.g.,  those  presented  57  mter  discharge  at  upstream 

in  the  Model  Rigor  section),  ice  momen-  md/  mid.reach/  and  dovmstream  end  0f 
turn  has  a  negligibly  small  effect.  SyS[em% 

For  example,  the  final  jam  thickness 
profile  in  Figure  42  is  within  10%  of  the 

equilibrium  thickness  value  for  that  water  flow  rate,  for  the  jam  with  an  initial 
thickness  of  1.3  m.  The  average  jam  thickness,  however,  is  1.5  m  or  2%  above  the 
equilibrium  value.  For  this  case,  a  steady-state  equilibrium  thickness  calculation 
would  have  proved  satisfactory.  When  unsteady  boundary  condition  data  are  not 
available,  it  may  be  necessary  to  use  a  simpler  steady-state  model  to  provide  jam 
thickness  estimates,  rather  than  using  the  fully  coupled  unsteady  model.  There¬ 
fore,  it  is  useful  to  identify  a  parameter,  or  set  of  parameters,  that  indicate  when  ice 
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Figure  57.  Water  discharge  at  upstream 
end,  mid-reach,  and  downstream  end  of 
the  system. 
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momentum  is  important  and  should  be  taken  into  account  for  determining  jam 
thickness  profile.  For  these  cases,  the  results  of  simpler  steady-state  models  should 
be  used  with  caution. 

The  force  balance  on  a  jam  at  static  equilibrium  can  be  represented  by  equations 
such  as  eq  19  or  150  with  the  nonuniform  terms  set  to  zero.  The  downstream-acting 
forces  of  the  water  shear  stress  on  the  jam  underside  and  the  component  of  gravity 
due  to  the  weight  of  the  jam  are  balanced  by  the  shear  resistance  at  the  banks 


2 

+  gS0Bs^  =  gs{  (l  -  p)(l  -  s;  )k0XKpr\2  .  (153) 

If  the  simplified  notation  of  a  +  b  =  c  is  used  to  describe  eq  153,  then  a/c  represents 
the  portion  of  the  jam  strength  (represented  by  the  bank  shear  resistance)  mobi¬ 
lized  by  water  shear  stress.  By  assuming  a  wide,  rectangular  channel,  the  first  term 
in  eq  153  can  also  be  written  as 


a  = 


fiU2B  jQS0gj2B 


=  /i 


8/o 


(154) 


Together  with  the  equilibrium  thickness  given  by  eq  25,  Table  2  was  generated  for 
jams  at  the  verge  of  stability  (equilibrium  thickness)  for  flow  in  channels  of  differ¬ 
ent  bed  slopes. 

Table  2  shows  that,  for  very  low  bed  slopes,  downstream  gravity  forces  are  mini¬ 
mal  with  the  water  shear  stress  mobilizing  most  of  the  jam's  strength.  For  the  smallest 
bed  slope,  when  the  discharge  was  increased  by  a  factor  of  2,  a/c  remained  very 
high  and  the  equilibrium  thickness  increased  by  24%.  At  the  highest  bed  slope 
shown,  water  shear  stress  only  engages  a  small  part  of  the  jam's  strength.  A  similar 
factor  of  2  increase  in  discharge  increases  a/c,  but  it  still  remains  low  and  only  a  6% 
change  results  in  the  equilibrium  thickness.  Smaller  bed  slopes  are  associated  with 
lower  water  velocities,  but  also  much  lower  jam  thicknesses.  As  bed  slope  increases, 
water  shear  stress  increases,  but  at  a  slower  rate  than  the  gravity  force  term  (b  in  eq 
153).  Once  a  jam  fails,  it  is  water  shear  stress  that  transports  the  ice  downstream. 
While  this  finding  is  not  entirely  intuitive,  it  means  that  ice  momentum  effects  be¬ 
come  more  important  as  the  a/c  ratio  increases,  i.e.,  for  smaller  bed  slopes  and 
smoother  channels. 

Consider  further  the  a  and  c  terms  and  their  ratio  in  eq  153  for  a  bed  slope  of 
0.0005  given  in  Table  2.  Equation  154  shows  that,  for  a  given  discharge,  the  a  term 
has  a  set  value.  Equation  153  balances  the  forces  by  adjusting  to  the  value  of  equi¬ 
librium  thickness,  giving  a  value  for  a/c  at  the  limit  of  stability.  But,  for  example,  if 
the  thickness  had  attained  a  value  greater  than  the  equilibrium  level  (say  1.60  m). 


Table  2.  Ice  parameters  for  channels  at  different  bed  slopes. 


Bed  slope ,  S0 

r\e(]  at  100  m3/s 

a/c  at  100  m3/s 

T\ec]  at  200  m3/s 

a/c  at  200  m3/s 

0.000005 

0.21m 

0.960 

0.26  m 

0.968 

0.00005 

0.49  m 

0.829 

0.60  m 

0.861 

0.0005 

1.47  m 

0.432 

1.69  m 

0.509 

0.001 

2.30  m 

0.277 

2.57  m 

0.352 

0.0025 

4.74  m 

0.121 

5.02  m 

0.171 

73 


eq  153  would  no  longer  be  in  balance  and  would  represent  a  condition  of  greater 
jam  stability.  In  this  case,  the  a/c  ratio  decreases  because  the  a  term  remains  con¬ 
stant.  As  water  discharge  increases,  the  jam  should  remain  stable  for  a  longer 
period  prior  to  shoving  and  thickening,  thus  resulting  in  lower  ice  velocities  and 
less  of  an  effect  of  ice  momentum  on  the  jam  thickness.  Also,  for  a  given  initial  jam 
thickness  (whether  at  the  limit  of  stability  or  not),  a  smaller  discharge  rise  has  less 
effect  on  jam  thickness  than  does  a  larger  one.  Finally,  a  larger  relative  discharge 
increase  has  more  effect  than  does  a  smaller  one.  For  instance,  ice  momentum  would 
be  expected  to  influence  the  final  jam  thickness  more  for  an  increase  from  100  to 
200  m3/s  than  an  increase  from  200  to  300  m3/ s. 

To  express  these  trends,  a  dimensionless  parameter  was  developed  that  includes 
initial  jam  conditions  (indicating  how  close  the  jam  is  to  the  limit  of  stability),  as 
well  as  the  relative  increase  in  discharge  expected.  This  number  is 

AQ"|  {  /i»2B _ 

Qin  J  (  8£S,  (1  -  p)(l  -  Si  )k0\KpT\2 

where  Qin  is  initial  water  discharge  and  AQ  is  the  expected  change  in  discharge. 
This  dimensionless  parameter  is  the  product  of  the  initial  state  of  stability  of  the 
jam  and  the  relative  discharge  increase  applied  to  cause  an  instability. 

Several  runs  were  made  with  the  fully  coupled  model  using  an  inflow  hydrograph 
that  rose  at  the  same  rate  as  the  baseline  inflow  hydrograph,  but  with  ten  different 
combinations  of  initial  discharge  and  discharge  increase  as  listed  in  Table  3. 

Runs  were  made  for  eight  different  bed  slopes  of  0.00005, 0.00008, 0.0001, 0.00025, 
0.0005,  0.00075,  0.001,  and  0.0025.  The  final  jam  thickness  profile  for  each  of  these 
runs  was  compared  to  the  equilibrium  jam  thickness  T|eq  for  the  final  discharge  as 
calculated  by  eq  25.  Average  values  of  jam  thickness  r\  were  calculated  for  each 
profile.  The  values  of  Q  were  plotted  against  T|  /  r|eq  in  Figure  58.  The  data  points 
delineate  a  line,  showing  combinations  of  channel  and  flow  conditions,  or  II  val¬ 
ues,  for  which  ice  momentum  significantly  affects  jam  thickness  (from  the  T|  /  T|eq 
value).  It  is  clear  that  ice  momentum  is  very  important  for  low  bed  slope  values, 
because  water  shear  stress  engages  the  greater  portion  of  the  jam  strength  for  these 
cases.  The  scatter  in  the  data  at  higher  values  of  T|  /  r|eq  is  most  likely  attributable 
to  the  highly  nonuniform  jam  thickness  profiles  for  those  cases.  Figure  59  shows  a 
plot  of  the  final  jam  thickness  profiles  for  cases  of  relatively  small,  medium,  and 
large  ice  momentum  effects  (Tl/^leq)-  A  smaller  ice  momentum  effect  results  in  a 
significantly  more  uniform  thickness  profile. 


Table  3.  Characteristics  of  various  inflow  hydrographs. 


Hydrograph  type 

Initial  Q  (m3/s) 

Final  Q  (m3/s) 

AQ  (m3/s) 

AQ/Q 

1 

100 

112.5 

12.5 

0.125 

2 

100 

125 

25 

0.25 

3 

100 

150 

50 

0.50 

4 

100 

175 

75 

0.75 

5 

100 

200 

100 

1.00 

6 

100 

250 

150 

1.50 

7 

100 

300 

200 

2.00 

8 

112.5 

200 

87.5 

0.78 

9 

125 

200 

75 

0.60 

10 

150 

250 

100 

0.67 
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Figure  58.  Dimensionless  jam  thickness  vs. 
ice  momentum  parameter  for  several  bed 
slopes. 


Figure  59.  Final  jam  thickness  profiles  for 
large ,  medium ,  and  small  ice  momentum 
effects. 


A  curve  such  as  the  one  presented  in  Figure  58  could  be  used  to  determine  the 
general  effects  of  ice  momentum  on  predicted  changes  in  water  discharge.  A  few 
simple  calculations  would  allow  a  determination  of  whether  a  static,  steady  flow 
model  is  satisfactory  or  if  a  fully  coupled  model  such  as  the  one  presented  here  is 
required  to  provide  accurate  jam  thickness,  and  thus  water  surface,  profiles. 

Effects  of  hydrograph  shape  on  jam  thickening 

The  experiments  discussed  in  the  previous  section  show  that  the  relative  state  of 
initial  stability  of  a  jam,  in  conjunction  with  the  relative  change  in  water  discharge 
expected,  determines  whether  ice  momentum  significantly  affects  jam  thickness 
profile.  Other  factors  may  also  determine  when  ice  momentum  is  significant.  Field 
observations  have  shown,  in  general  terms,  that  fast  rises  in  discharge  because  of 
heavy  rain  on  frozen  ground  often  result  in  severe  ice  jams  producing  devastating 
flooding.  Commensurably,  gradual  warming  trends  with  slower  discharge  increases 
and  higher  air  temperatures,  which  cause  the  jam  strength  to  weaken,  often  result 
in  less  dynamic  or  damaging  jams.  The  experiments  discussed  in  the  Laboratory 
Experiments  section  attempted  to  simulate  short  duration  peak  flow  events  by  in¬ 
creasing  the  water  discharge  and  then  decreasing  it  again  after  a  short  period.  The 
effect  of  these  short  duration  peak  flows  was  to  stall  the  progression  of  the  shoving 
and  thickening,  resulting  in  only  partially  thickened  jams.  For  those  cases,  the  ice 
momentum  may  be  important,  but  the  peak  discharge  may  be  attenuated  greatly 
by  the  time  the  disturbance  travels  the  length  of  the  entire  system. 

A  series  of  experiments  was  conducted  to  assess  the  effects  of  the  shape  of  the 
inflow  hydrograph  on  the  resulting  jam  thickness  profile  and  to  establish  whether 
ice  momentum  continues  to  be  a  consideration  in  these  cases.  Runs  were  made 
with  various  rates  of  rise  on  the  rising  limb  of  the  hydrograph.  The  baseline  runs 
had  a  rate  of  rise  of  2.5  m3/s  for  each  time  step  of  0.5  minutes.  For  these  runs,  the 
time  step  was  held  constant,  but  the  boundary  condition  file  describing  the  inflow 
hydrograph  was  modified  for  rise  rates  of  between  10  and  0.3125  m3/s  for  the  same 
0.5-minute  time  step.  This  modification  resulted  in  times  to  peak  £p  (time  for  the 
hydrograph  to  increase  from  100  to  200  m3/s)  of  5, 10,  20, 40,  80,  and  160  minutes. 
The  inflow  hydrographs  are  shown  in  Figure  60.  Time  to  peak  £p  was  related  to  the 
base  time  defined  as  the  time  taken  for  a  disturbance  to  travel  the  length  of  the 
simulated  channel.  From  eq  152 
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where  LT  is  the  total  length  of  the  mod¬ 
eled  system.  Runs  were  made  for  the 
values  of  fp  given  above  and  three  dif¬ 
ferent  bed  slopes — 0.00005,  0.0005,  and 
0.001.  A  functional  relationship  between 
the  ratio  of  average  final  jam  thickness 
to  equilibrium  thickness  T|  /  T|eq  and 
fp/  fb  is  evident  in  Figure  61.  As  the  time 
to  peak  fp  increases,  the  average  final  jam 
thickness  r[  approaches  that  predicted 
by  equilibrium  theory  Tjecj.  The  effects  of 
fp  on  the  shape  of  the  final  jam  thickness 
profile  is  demonstrated  in  Figure  62, 
which  shows  the  final  profiles  for  times 
to  peak  of  5,  20  (standard),  and  80  min¬ 
utes  for  the  bed  slope  of  0.0005.  As  fp 
decreases,  the  changes  in  water  dis¬ 
charge  for  each  time  step  increase,  result¬ 
ing  in  larger  gradients  of  the  variables 
in  the  upstream  reaches.  The  gradient 
terms  (e.g.,  dr\/dx  )  are  not  included  in 
the  equilibrium  thickness  formulation  of 
eq  25  and  play  an  important  part  in  the 
full  ice  momentum  equation.  Also,  the 
large  changes  in  water  discharge  easily 
overpower  the  stability  of  the  jam,  result¬ 
ing  in  significant  ice  velocities  and  ice 
momentum  effects. 

Experiments  were  also  conducted  to 
investigate  the  effects  of  the  sustained 
time  at  the  peak  flow  fs.  For  these  runs, 
the  standard  inflow  hydrograph  rate  of 
rise  of  2.5  m3/s  increase  per  time  step 
was  used,  but  the  time  at  which  the  wa¬ 
ter  flow  was  held  at  200  m3/s,  before  be¬ 
ing  reduced  back  to  the  initial  value  of 
100  m3/s,  was  varied.  The  sustained 
times  simulated  were  0  (instantaneous 
peak),  10,  20,  40,  80,  and  °°  minutes  (no 
discharge  decrease)  as  shown  in  the  in¬ 
flow  hydrographs  depicted  in  Figure  63. 
These  runs  were  conducted  at  the  same 
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Figure  60.  Inflow  hydrographs  for  various 
times  to  peak. 
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Figure  61.  Dimensionless  jam  thickness  vs. 
tp/t b  at  three  bed  slopes. 
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Figure  62.  Effect  of  tp  on  final  jam  thick¬ 
ness  profile  shape. 


three  bed  slopes  as  above. 

Figure  64  shows  relationships  between  the  ratio  of  final  average  thickness  to 
equilibrium  thickness  and  tj  fb.  For  the  case  of  the  infinite  sustained  time,  a  nomi¬ 
nal  value  of  fs/fb  (beyond  which  there  were  no  changes  in  jam  thickness)  was  cho¬ 
sen  for  plotting  purposes.  Again,  the  effect  of  bed  slope  is  clearly  seen.  It  is  signifi¬ 
cant  that,  as  the  value  of  ts/ fb  reaches  1,  the  effects  of  further  increases  in  the 
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Figure  63.  Inflow  hydrographs  for  various 
times  of  sustained  flow. 
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sustained  time  on  the  final  thickness  pro¬ 
file  become  negligible.  The  plot  does 
show,  though,  that  even  with  an  instan¬ 
taneous  discharge  peak,  the  effects  of  ice 
momentum  are  important. 

Figure  65  presents  the  changes  in  the 
final  thickness  profile  with  increased  ts 
for  the  0.0005  bed  slope.  The  increase  in 
ts  results  in  a  slight  smoothing  and  fur¬ 
ther  thickening  of  the  downstream 
reaches.  Even  for  the  instantaneous 
peak,  however,  the  jam  thickness  in  the 
upstream  reaches  is  significantly  greater 
than  the  equilibrium  value  of  1.7  m. 
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Figure  64.  Dimensionless  jam  thickness  vs. 
tftfr  at  three  bed  slopes. 


SUMMARY 

It  is  well  known  that  ice  jams  are  in¬ 
herently  unsteady  events,  in  which  mov¬ 
ing  ice  is  brought  to  rest  as  accumula¬ 
tions  that  shove  and  thicken  in  accor¬ 
dance  with  changing  forces  exerted  by 
water  flow,  accumulation  weight,  and 
bank  roughness.  These  processes  are 
even  more  unsteady  when  a  jam  col¬ 
lapses,  plows  downstream,  and  possibly 
reforms.  Not  well  known,  at  least  in  any 
quantitative  way,  is  how  the  inherent 
unsteadiness  of  those  processes  affects 
the  jam  formation  and  thickness  profile. 
In  particular,  virtually  nothing  is  known 
about  the  effect  of  ice  momentum  on  jam 
formation  and  thickening.  This  study 
presents  the  first  formulation  and  exami¬ 
nation  of  the  fully  coupled  dynamic  na¬ 
ture  of  the  unsteady  processes  associated 
with  jam  formation.  It  does  so  by  means 
of  numerical  simulations  and  laboratory 
flume  experiments.  The  experiments 
also  show  that  equilibrium  thickness  for¬ 
mulations  consistently  underestimated 
measured  jam  thickness. 

The  numerical  simulation  uses  the  full  one-dimensional  unsteady  flow  equa¬ 
tions  (conservation  of  mass  and  momentum)  for  water  flow,  ice  movement,  and 
jam  formation.  A  unique  aspect  of  the  simulation  is  that  the  equations  are  solved  in 
a  fully  coupled  manner.  The  simulation  model  is  shown  to  be  robust,  versatile,  and 
accurate  in  calculating  jam  thickness  profile  under  a  variety  of  initial  jam,  flow,  and 
channel  conditions.  The  effects  of  the  shape  of  the  inflow  water  hydrograph  on  jam 
thickness  profile  are  investigated  and  shown  to  be  important.  A  dimensionless  pa- 


Figure  65.  Effect  ofts  on  final  jam  thickness 
profile  shape. 
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rameter  is  used  to  delineate  general  conditions  for  which  ice  momentum  signifi¬ 
cantly  affects  jam  thickness  profile.  The  parameter  takes  into  account  initial  ice 
conditions  and  the  expected  flow  changes. 


CONCLUSIONS 

The  study  led  to  the  following  principal  conclusions: 

1.  The  flume  experiments,  which  provide  a  detailed  description  of  the  observed 
processes  of  jam  formation,  failure,  and  evolution,  brought  to  light  two  jam  failure 
and  reformation  mechanisms — progressive  and  complete.  Progressive  jam  failure 
and  reformation  happens  with  lower  initial  discharge  and  lower  discharge  increases 
relative  to  the  discharge  needed  to  completely  destabilize  the  jam.  It  is  character¬ 
ized  by  a  smooth  progression  of  a  shoving  front  downstream  through  the  jam.  Com¬ 
plete  jam  failure,  then  jam  reformation,  occurred  for  initial  discharges  close  to  the 
discharge  necessary  to  completely  destabilize  the  jam.  It  is  characterized  by  the 
entire  jam  mobilizing  en  masse,  moving  downstream,  slamming  into  a  downstream 
barrier  (e.g.,  an  existing,  stationary  ice  cover,  or  an  ice  boom),  and  reforming. 

2.  The  flume  experiments  revealed  that  progressive  and  complete  modes  of  jam 
failure  and  reformation  result  in  measured  jam  thicknesses  that  exceeded  those 
predicted  using  prior  jam  formulations,  based  on  analyses  of  stationary  jams.  This 
finding  confirms  that  the  momentum  of  the  moving  ice  arrested  during  jam  forma¬ 
tion  produces  an  important  force  that  should  be  taken  into  account  when  estimat¬ 
ing  jam  thickness  for  many  conditions.  Prior  formulations  underestimate  jam  thick¬ 
ness  because  they  do  not  include  this  force,  or  the  interaction  of  the  water  and  ice. 

3.  The  numerical  simulation  model  provides  further  quantitative  information 
illuminating  the  effects  of  ice  momentum  on  jam  thickness.  Not  only  does  ice 
momentum  result  in  greater  average  jam  thicknesses  than  predicted  using  the  sta¬ 
tionary  jam  theory  of  prior  formulations,  it  also  induces  a  high  degree  of 
nonuniformity  in  a  jam  thickness  profile.  The  interaction  of  the  ice  movement  and 
water  flow  result  in  unsteady  variations  in  water  depth  and  velocity,  as  well  as  ice 
thickness  and  velocity,  throughout  the  entire  simulated  flow  and  jam.  It  is  these 
interactive  effects  that  result  in  nonuniform  jam  thickness  profiles.  For  this  case  of 
zero  ice  velocity  at  the  downstream  boundary,  the  profiles  are  characterized  by 
greater  thickness  in  the  upstream  reaches  where  ice  velocity  and,  thus,  momentum 
are  greatest.  The  use  of  a  fully  coupled  solution  technique  preserves  this  interac¬ 
tion  of  the  variables.  The  loosely  coupled  solution  results  in  some  averaging  or 
smoothing  of  the  variables  and  their  effects  upon  each  other.  As  the  time  step 
decreases  or  the  number  of  ice  calculation  cycles  increases,  the  results  of  the  loosely 
coupled  model  approach  those  of  the  fully  coupled  model. 

4.  The  dimensionless  parameter  Q  =  (a/c)(AQ/Qin)  is  useful  for  delineating  con¬ 
ditions  when  ice  momentum  should  be  taken  into  account  for  jam  thickness  esti¬ 
mation.  In  £1,  a  =  (fiu2B)/S/  and  c  =  ^si(l-p)(l-si)/c0^KpT|2,  so  that  (a/c)  represents  the 
portion  of  the  initial  jam  strength  mobilized  by  the  water  shear  stress  on  its  under¬ 
side.  The  parameter  relates  the  ratio  of  average  jam  thickness,  as  determined  by  the 
fully  coupled  model,  to  steady  equilibrium  thickness  determined  from  stationary 
jam  models.  The  relationship  is  useful  for  establishing  when  changes  in  flow  con¬ 
ditions  (i.e.,  hydrograph  properties)  will  destabilize  a  jam  and,  through  ice  impact, 
affect  jam  thickness  profile.  The  parameter  delineates  the  conditions  when  a  fully 
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coupled  model  of  jam  evolution  is  needed  to  accurately  predict  jam  thickness. 

5.  The  shape  of  an  inflow  water  hydrograph  and  its  peak  flow  are  important  in 
estimating  jam  thickness  because  they  determine  whether  jams  fail  and  reform  pro¬ 
gressively  or  completely.  Fast-rising  hydrographs  with  attendant  rising  water  lev¬ 
els  result  in  larger  gradient  terms  across  a  computational  reach.  The  gradient  terms 
(including  ice  momentum)  result  in  greater  predicted  jam  thickness.  They  are  not 
present  in  equilibrium  thickness  formulations,  which  thereby  underpredict  thick¬ 
ness  and  water  levels.  The  time  that  a  flow  remains  at  its  peak  value  is  also  impor¬ 
tant.  Instantaneous  peaks  are  attenuated  as  they  travel  downstream,  resulting  in 
lower  local  discharges  and  lower  stresses  on  the  underside  of  the  jam,  leading  to 
thinner  jams. 


RECOMMENDATIONS  FOR  FUTURE  RESEARCH 

The  flume  experiments  and  numerical  model  bring  to  light  several  areas  where 
deficiencies  in  knowledge  and  formulation  persist. 

Several  parameters  are  necessary  to  characterize  the  properties  of  ice  accumu¬ 
lated  in  a  jam.  Among  the  properties  are  angle  of  internal  resistance  <j>,  coefficient  of 
friction  of  ice  on  ice  X,  and  the  lateral  expansion  coefficient  Practical  measure¬ 
ments  of  <|)  have  been  obtained  for  particulate  materials,  with  an  estimate  given  as 
the  dry  angle  of  repose  of  the  material.  No  known  experiments  have  been  con¬ 
ducted,  however,  to  determine  suitable  values  of  X  or  for  ice.  Further  laboratory 
experiments  to  define  them  are  necessary. 

The  full  effects  of  ice  momentum  have  not  yet  been  fully  identified  because  of 
the  influence  of  ice  velocity  on  the  parameters  mentioned  above.  Of  particular 
importance  is  the  temporal  variation  of  these  parameters  as  an  ice  jam  fails  and 
mobilizes,  during  which  the  forces  on  a  jam  continually  change.  The  value  of  is 
related  to  <[>  and  has  been  shown  to  be  adequately  described  by  Mohr-Coulomb 
theory  for  a  stationary  jam.  The  Rankine  states  of  active  and  passive  pressure 
described  using  Mohr-Coulomb  theory  are  intended  as  point  values  at  the  instant 
of  particulate  material  failure.  Questions  surround  the  value  of  the  passive  pres¬ 
sure  coefficient  once  a  jam  fails  or  is  moving  downstream.  A  similar  argument  goes 
for  the  value  of  coefficient  of  friction  of  the  ice  along  the  shear  boundary  at  the 
banks  X.  Certainly,  if  this  coefficient  is  developed  in  a  way  similar  to  that  for  the 
simple  static  coefficient  of  friction,  there  must  be  changes  that  occur  as  the  material 
comes  into  motion.  It  is  likely  that  the  coefficient  would  reduce  once  a  jam  moves, 
thereby  resulting  in  a  lower  resistance  to  downstream  movement.  Thicker  jams 
would  result  when  the  moving  ice  is  finally  arrested  because  bank  friction  is 
reduced. 

Several  modifications  to  the  model  are  possible  to  make  the  simulations  more 
realistic.  The  modifications  include  the  ability  to  use  actual  cross-section  geometry, 
and  changes  in  roughness  coefficients,  channel  width,  and  ice  parameters  with  dis¬ 
tance.  Extending  the  model  to  two-dimensional  coordinates  would  not  only 
increase  the  computational  time,  but  would  also  introduce  difficulties  in  defining 
the  failure  mechanism.  The  model  developed  in  this  study  assumes  that  ice  param¬ 
eters  are  satisfactorily  given  as  bulk  properties.  A  more  suitable  advance  might  be 
the  adaptation  of  the  model  to  simulate  branched  systems,  where  jam  formation  in 
one  channel  would  affect  water  and  ice  movement  in  another. 
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The  very  local  problem  of  jam  failure  also  is  a  candidate  for  future  research.  The 
vertical  stress  level  throughout  the  thickness  of  a  jam  varies  linearly  when  the  jam 
is  at  rest.  This  stress  description  is  the  basis  for  the  determination  of  K p  and  other 
jam  parameters.  When  water  shear  stress  exists  beneath  the  jam,  however,  the  stress 
levels  within  the  jam  change.  No  information  exists  about  these  stress  levels  and 
whether  they  have  an  effect  in  setting  up  the  initial  instability  within  the  jam.  The 
vertical  description  of  jam  failure,  whether  the  jam  initially  fails  at  the  bottom  or  at 
the  water-surface  level,  where  stress  is  assumed  to  be  the  maximum,  does  not  exist. 

The  highly  dynamic  nature  of  ice  jamming  precludes  direct  measurement  of  thick¬ 
ness  and  ice  velocity  during  formation,  evolution,  and  failure.  Equally  difficult  to 
obtain  are  time-histories  of  water  discharge  and  depth  at  several  locations  within  a 
reach  experiencing  jamming.  A  concerted  effort  needs  to  be  undertaken  to  obtain  a 
complete  set  of  field  data  for  an  ice  jamming  event,  including  water  depths,  water 
velocities,  ice  thicknesses,  and  ice  velocities  at  several  locations.  Another  avenue  of 
model  verification  and  validation  would  be  a  detailed  physical  model  study  of  ice 
jamming,  going  beyond  the  level  of  the  flume  experiments  described  in  the  Labora¬ 
tory  Experiments  section.  Such  a  data  set  will  allow  the  verification  and  comparison 
of  unsteady  ice  models,  such  as  the  one  presented  in  this  study. 
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APPENDIX  A:  DISCRETIZATION  OF  FUNCTIONS  AND  DERIVATIVES 


The  Newton-Raphson  iteration  procedure  required  not  only  the  discretization 
of  the  equations  of  motion,  but  also  the  partial  derivatives  of  these  discretizations 
with  respect  to  the  solution  variables.  For  example,  the  conservation  of  water  mass 
equation  is  discretized  as 
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The  conservation  of  ice  mass  equation  F2  is  given  as 
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(All) 


and  the  partial  derivatives  are  evaluated  as 


=  ^ _ *#1  t  6j 

9rij  2A  t  Ax  2 


_’>h 


c  2  Ax 


c'  =  ^  =  0 

OMj 


d'  =  ¥*=?i 

dVj  2 


L  ei  (rif++il  ~~T>r+1)+(i~ei  ^r+i  - 

>  Ax  Ax 


S.&--0 


3F2  _  1  |  xBj  |  6,  9i  K+l1  -  •»)" +1 )  +  (!  -  ei  )(^f+i  -  ) 

"V~.  T  Al  A _  O  A 


^rjj+1  2A  t  Ax  2 

s'=^-=o 


^/=_a^_  =  eL  QiK-tf  -<+1)+(i~9i)K+i -d" )  +^eL 

3\)j+1  2  Ax  Ax 


The  conservation  of  water  momentum  equation  is  discretized  as 
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The  conservation  of  ice  momentum  equation  is  discretized  as 
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and  the  partial  derivatives  are  evaluated  as 
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The  reach  averages  of  the  dependent  variables  and  their  combinations  are  given  as 
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of  their  effects  on  jam  thickness.  Ice  momentum  should  be  taken  into  account  for  most  jams  because  it  leads  to 
significantly  thicker  jams  and  affects  the  thickness  profile.  A  useful  dimensionless  parameter  is  identified  for 
generalizing  this  finding. 
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